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ABSTRACT 


An equation solving package based on the skyline technique is 
developed for substructure analysis, as well as an assembly and a co- 
ordinate transformation scheme for the same purpose. Two plane frame 
substructure analysis programs, SISAPF and MUSAPF, are developed. The 
former is based on a single-level substructure scheme, and the latter is 
based on a multi-level substructure scheme. 

A sample structure is partitioned in several ways. The core 
space parameters, and the CPU time requirements for the various stages 
of solution are discussed. A large saving in core space requirements, 


and a significant saving in CPU time are achieved. 
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CHAPTER 1 - INTRODUCTION 


1.1. Introductory Remarks 


Matrix methods of structural analysis are based on representing 
the actual structure by discrete structural elements interconnected at a 
set of nodes. The matrices representing the elastic properties of the 
individual elements are assembled to model the total structure. This 
mathematical model is used to determine the forces and displacements at 
the joints. When these are known, the conditions within each element 
may be determined. 

The substructure method of analysis is a natural development 
of this concept. A structure can be partitioned into a number of sub- 
structures, which are considered as complex structural elements. These 
substructures are assembled along the inter-boundary nodes (Fig. 1.1). 
Once the inter-boundary system is solved, the values of the unknowns 
associated with the internal degrees of freedom of the complex struc- 
tural elements can be obtained. 

The concept can be extended to model what might be called a 
hierarchy of substructures [5,4]. This means that a substructure may 
itself be an assemblage of small substructures, which may be assemblies 
of yet smaller units, until one reaches simple structural elements (Fig. 
Wey: 

There are two common methods of substructure analysis. The 
first method [11] obtains the boundary matrices through full release of 
the internal degrees of freedom. The second method [5,12] achieves the 
same result through decomposition, sometimes called partial release of 


the internal degrees of freedom. These two approaches are discussed 
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Inter-boundary Nodes 


Internal Node 


Simple Structural Element 


Figure 1.1 Substructures 


Figure 1.2 A Hierarchy of Substructures 
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in the following section. 


1.2 Full Release vs Partial Release 
Consider Eq. 1.1 which represents the load displacement relation 


of a substructure, without imposed boundary displacements, 
ERP ey = eR} (171) 


where {r} is the vector of nodal displacements, {R} is the vector 
of nodal forces, and [K] is the substructure stiffness matrix. Assume 
that the degrees of freedom can be partitioned into internal degrees of 
freedom, denoted by subscript i, and inter-boundary degrees of freedom, 


denoted by subscript b. Eq. 1.1 can be put in the form 


Kiy Kb R, 
= (1.2) 
Ka: K ry Re 
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and 


* = 
{Ri Je = {Rpt - [Kid [K. 5] : {R.} (1.3d) 


If boundary conditions are imposed on rare the vector tr} can be 
determined from Eq. 1.3b. Substituting this vector into Eq. 1.3a, the 
vector {r.} can also be obtained. This is calied the ‘full release' 
approach [11,12]. 

Another approach is to consider the component matrices of Eq. 
1.2 as matrix elements, and carry out a simple Gaussian elimination on 


the first row, to get 


x x 
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The inter-boundary partitions of Eq. 1.4 are identical to those of Eqs. 
1.3c, and 1.3d. Eq. 1.4 is formed by a decomposition of Eq. 1.2 up to 
the last row of the internal degrees of freedom [12,15]. This latter 
approach is the ‘partial release' method, sometimes cailed the "Aitkin" 
decomposition [5}. 

Williams [15] has evaluated the number of numericai operations 
required to obtain the inter-boundary partitions using either approach. 
Table 1.1 simplifies nis results for one Baas OG loading. It can be 
seen that the partial release approach involves fewer numerical opera- 


tions. Moreover, this table is very conservative in the case of partial 
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release since it assumes fully populated matrices. Because of the 
greater efficiency of the partial release technique the full release 


technique is seldom used and will not be discussed further herein. 


1.3 Purpose and Scope 


The objectives of this thesis are: 
1) To develop an efficient equation solver to deal with the various 
stages of substructure analysis. 
2) To develop a flexible single-level substructure analysis program, 
and a practical highly flexible multi-level substructuring scheme. 
3) To explore the application of both single-level substructuring, and 


multi-level substructuring to planar framed structures. 


1.4 Qutline of Contents 

Chapter 2 contains the theoretical background for, and the 
development of, an equation solving package suitable for substructure 
analysis. General techniques of equation solving are outlined. The 
"skyline' method is developed for substructures, and the required 
algorithms are derived. The efficiency of the ‘skyline’ method is then 
compared to other methods in a qualitive way to indicate its advantages. 

Considerations of assembly on coordinate transformation of 
the inter-boundary partitions are dealt with in Chapter 3. The logic 
flow and concepts of the single-level substructure analysis program 
SISAPF, followed by the conceptual development and the logic flow of the 
multi-level substructuring program MUSAPF, form the contents of Chapter 4. 


In Chapter 5, an example problem is handled in several ways, 
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and the effects of the number, and the number of levels, of substructures 
on execution time are discussed. In Chapter 6 conclusions are drawn, 


and possible future developments are outlined. 
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CHAPTER 2 - SOLUTION OF EQUATIONS 


2.1 Introduction to Solution of Equations 


The equation solver is an essential part of any structural 
matrix analysis program. Several techniques have evolved since the 
advent of digital computers; all of them aimed at efficiency, since the 
CPU time consumed in equation solving usually constitutes 20% to 50% of 
the total time required to process a problem [9]. Reference 9 presents 
a review of these methods, and Reference 13 contains valuable informa- 
tion on the concepts involved in decomposition techniques. Some of the 
basic methods are reviewed in the following sections (Sect. 2.2 to 2.7). 
The skyline algorithms, developed by the author for implementation in 
the substructure analysis programs, form the subject material for the 


remainder of the Chapter. 


2.2 Standard Gaussian Elimination 


Consider the system of equations 
bad Ou t=) {1b} (2a) 
in which [a] is a coefficient matrix, {x} is the vector of unknowns, 


and {b} is the (known), so-called, right hand side. The algorithm for 


elimination, with equation i as the pivotal equation, will yield 
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The nth elimination will yield 


on F by fain | (2.3) 


The unknowns Xap? Xng? Xz> are found by backsubstitution, using the 


equation 
n 
me =O(pumen tS’ 927 x2 )}Vaw levi =tnSsl oe) (2.4) 
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Introducing symmetry, which is a property of most structural 
matrices, aiy = ayy? and Eqs. 2.2a, and 2.2b become 

eke ad ak ie eo 

Siren tn (2.5a) 
bate Dust-vae. beyae. 6) =o 0 Fe A (2.5b) 


Further, it should be noted that structural matrices for stable struc- 
tures are always positive definite and well posed, hence pivoting is 


rarely justified [9] and will not be discussed. 


2.3 Gaussian Elimination as Matrix Decomposition 


Gaussian elimination can be interpreted as a decomposition of 
matrix [a] into upper and lower triangular matrices [a],; and [a]\> such 


that 
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The solution of Eq. 2.1 may be represented as 


OG eek be Gier (2.7a) 


{X} fale {Y} (227) 
The coefficients of [a], are those obtained from Eq. 2.2a. 
The diagonal components of [a], are unity. The evaluation of {Y} 
corresponds to the evaluation of {b} by Eq. 2.2b or 2.5b, and the 
evaluation of [a], {Y} is the backsubstitution process described by 
Eq. 2.4. The Gaussian elimination process is thus a decomposition which 
proceeds by rows. It is possible to carry out the decomposition by 
columns, which is called Cholesky decomposition. For a symmetric matrix 


either decomposition can be put in the form 
a zs 
[a], [0] [a], = [al (2.8) 


where the diagonal elements of [a], in Eq. 2.8 are equal to unity, and 


[D] is a diagonal matrix. This will be discussed in detail subsequently. 


2.4 Banded Algorithms 


It is possible in an unbranched problem to number the nodes, 
such that the non-zero components of the coefficient matrix are clustered 


about the main diagonal. This type of matrix is called "banded', and 
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Figure 2.1 Symmetric Banded Matrix 


(a) Gauss (b) Cholesky 


Figure 2.2 Progress of Decomposition [13] 
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the non-zero components of the matrix may be represented schematically 
as shown in Fig. 2.la. ‘The limits in Eq. 2.5a become k = j, j tm- 1, 
and J = i+ 1, 1 +m- 1 in which m is called the ‘half band width! and 
includes the main diagonal. Further, since the matrix is symmetric, 
only the half band need be stored in a rectangular matrix, see Fig. 
2.1lb. Eqs. 2.5, adjusted for this type of storage will become 
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Fig. 2.2a shows the progress of a banded Gaussian elimination algorithm, 
where the shaded area is the space affected by elimination with equation 
i. Fig. 2.2b shows the progress of a banded Cholesky algorithm, and the 


shaded area is the space affected by decomposition of column i [7]. 


2.5 Frontal Methods 

These methods are Gaussian elimination schemes in which a row 
is eliminated as soon as all its components are assembled. The order of 
assembly of the elements then devennsetne maximum band width arising in 
the solution algorithm. This method, despite its elegance, requires a 
complicated housekeeping system. The technique is effective in large 
systems, where efficient nodal numbering poses difficulties, and where 
the perepherial processing of the problem is more time consuming than 
solving the equations [9]. However, the maximum band width required in 


a wave front solution is the same as that required with an efficient 
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band solution. 


2.6 Sparse Matrices 


In sparse matrices, where the band is extremely irregular, 
pivotal rows will contain large percentages of zero components. These 
components are inactive (Fig. 2.3a), and can be skipped during numerical 
Operations. However, in the case of a banded algorithm, they must be 
stored, unless the programmer introduces a separate addressing array. 

In a column, zero components below the first non-zero component may in 
general be expected to assume non-zero values during decomposition (Fig. 
2.3a). A large saving in storage requirements can be achieved if only 
the active columns of a matrix are stored. The active column associated 
with degree of freedom j is defined, herein, as that portion of column j 
that begins at the first non-zero component and ends at the diagonal 
Gridericesb): 

Columnwise decomposition is advantageous for the type of 
storage described above. In a substructured problem, where partial 
decomposition is required, the inter-boundary nodes should be grouped at 
one end of the nodal numbers, as indicated in Eq. 1.2, to avoid inter- 
changing rows and columns before decomposition and after backsubstitution. 
Optimal nodal numbering, which gives a minimum band width is no longer 
feasible, and the resulting substructure stiffness matrix is sparse, as 
illustrated in Fig. 2.4. 

For the reasons stated above, it was decided to use a column- 
wise decomposition scheme for the substructure analysis programs SISAPF, 
and MUSAPF. The method used is often referred to as the ‘skyline method’ 


[3,2], or a variable band width Cholesky method [13]. 


sor vueetazdus ‘ni PSP) hing yontts tie Shinty We Rt 


i Ee F aL ~ ‘a f ah i 
y LS, ' e ; i > i a ry 


setopore! ae ay beta any soul ie 


y alin 

§ rant” CTO Og RD on ~<apomnosngy AIRE anv 
io Pee aye uty eae uF was bap (ee Bs ie a ‘se $1 
Fea yods: audits ett panoad reed 1 ag 83 i at 
(ferret (tha seibbe tah hage B 2B aiiOren, rei ngid 
HL (RO SHS, Cen IY Ay ai sted cama 
als OT “Dogars wl “i 2oud ey via Agi omyaile- 3 
ny) Veatdgiesa dake eapehbidupas beat? fag 
| Bena oe mae: be Mi an ome 


; 
‘ H Ws rs 
L mules Ao-Anp 08 A iA us yee esa Pins Part > ve ho 
hs 
HINA Th “aris Ps abies Fie sea (salon. igiiey! vas. 


Pi i eu : 
| ¥ ¥, bit +h 


FAs rei, ie va 


asaya e ale Bey Las lem 
(erihia so amore y my ava she 
! A baquiwn oh leads Sahan: Cpl event: ‘pay Ost pon hun 


agent ‘bine 'od Fo Be ak modi te. seitiuns anuiiry anid ite 


fe PA SE Ura Os eer eae ond 


vengiter| oh wf tude "pew maak gibi vai one” Dvd ence tbian = 

7) “at wnged: a4 -9T4 aoa EL sui inane pnit tien gis, tins ‘ih bing ; 
oh Ma Mj U8 Loma bsraatewtt 1 
en ae v! habt ib aie 2, Svode Sainte. 2 erae net sett yo ane) I ‘a 
e TASS qiexgy 34 oh ay hging’ “wd wad, alta’ “Or. apse inet be rr ai : a 
‘heititwe ant Aine ont 2000) Lite S15 nine heztboltan: alt 4a bow r. " 


es 


“a LEC), borat leat pote D6 4106 tri (6 Ae 1 oS) 


eae —. . ; i é nae 
f > ‘ ® 


Active Column 


4 
| Active Row 


(a) (b) 


Figure 2.3. Advantage of Columnwise Storage 
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Figure 2.4 Sparseness Due to Substructure Nodal Numbering 
(Nodes 7, 8, 9, and 10 are the inter-boundary nodes.) 
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2.7 The Skyline Method 
This method has been described by Wilson and Bathe [2], and by 


Felippa [3]. It has been used in some large scale finite element pro- 
grams such as NONSAP [1], where it has been implemented for out of core 
solution of equations. It has two main features; the coefficient matrix 
storage scheme, and the columnwise decomposition algorithm. These are 


discussed below. 


2.7.1 The Coefficient Matrix Storage Scheme 


Since a structural stiffness matrix is symmetrical, only the 
upper triangle, including the main diagonal need be stored, and then 
only the active columns need be stored. The active column height for 
any degree of freedom is the difference between this degree of freedom 
and the lowest degree of freedom that appears in any element associated 
with it, see Fig. 2.5. It is tacitly assumed that all degrees of freedom 
associated with any one element are inter-active. By looping over all 
elements, determining the difference described above, and updating 
column heights as may be necessary, the skyline of the structure stiff- 
ness matrix is established. The skyline is defined as a bound on the 
active column heights (Fig. 2.6a). 

The active columns are oer in a one-dimensional array, see 
Fig. 2.6b. Access to any component is ensured if the addresses of the 
diagonal components are known. These are obtained by adding up the 
column heights, and storing the addresses of the diagonal components in 
a separate array, MAXA (Fig. 2.6c). For example, access to component 
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Figure 2.5 Column Heights for An Element (HT79 = Column height of degree 
of freedom 10 due to connectivity of element shown in (a)) 
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Figure 2.6 Storage Scheme 
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and A is the one dimensional storage array. 


2.7.2 Columnwise 'Cholesky' Decomposition 


This summary of the columnwise decomposition has been largely 
drawn from Reference 2. In the following the stiffness matrix will be 
denoted by [K], the load vector by {R} , and the displacement vector by 

{U} . The matrix corresponding to [a]! of Eq. 2.8 will be denoted by 
[L]. The stiffness matrix is symmetric, and if enough restraints are 
imposed on the system, it is positive definite. Hence, Eq. 2.1 can be 


put in the form 


[kL] [0] (L]' tuk = {Rr} (2.11) 
or 

PETeLG] StU Ser} (2.12a) 
where 

fe] .=,0pg-te1! (2.12b) 


Matrices [L], [G], and [D] are not explicitly obtained. Instead, the 
coefficients of [K] are modified. Eqs. 2.12a, and 2.12b yield the 
following algorithm, where the subscripts apply to a square matrix, and 
m. is the row subscript of the first non-zero element appearing in 
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Matrices [G], fae and [D] are formed by Eqs. 2.13a, 2.13b, and 2.13c 


respectively. The solution is then obtained as 


(0 aot eae: (2.14a) 
where 


{V} 


ieee re Ri (2.14b) 
Eq. 2.14b yields the forward substitution algorithm described by 
ey SIN ae Leet er ji =es 0 (2.15a) 
In coding this algorithm, the vector {V} overwrites {R} resulting in 
Roa Reg ee Kae Res Jo=ee5010 (2555) 


Eq. 2.14a describes the backsubstitution process and can be put in the 
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where the superscript indicates that variable is computed in backsub- 
stitution step number (n - superscript), and the final values of V of 


Eq. 2.17 are the vector {U} of Eq. 2.14a. 


2.8 Symbolic Development of Cholesky Decomposition for Substructure 
Analysis 
In order to determine the effective inter-boundary stiffness 
matrix [kK] by the partial release method it is required to reduce Eq. 


1.2 to the form 


* — kk 
a 0K U, R. 
: (2.18) 
* * 
Pd U, Re 


where [6] is an upper triangular matrix, and i and b denote the internal 
and inter-boundary degrees of freedom respectively. It will be shown 
that [K]” of Eq. 2.18 is identical to [K] of Eq. 1.3c. Since [K] is 
symmetric it can be decomposed to the form [L~] [D7] eye in which 

[D~] is composed of 4 submatrices, where the two off-diagonal submatrices 


are null matrices, and one of the partitions must have enough restraints 
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to stabilize the system. To demonstrate this, assume that the matrix 


product [L~] [D7] rosy! can be written in the form [12] 
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al pe One 
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Then the product of the matrices in Eq. 2.19 may be written as 
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If the matrix of Eq. 2.20 represents the stiffness matrix, the submatrices 


must identify with those of Eq. 1.2. This establishes the requirements 
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which is of the required form. Eq. 2.2la indicates that Toe and [D] 
may be determined from a standard decomposition of [K..]. Once this has 


been accomp1l ished iu may be determined from Eq. 2.21b as 
Ghd” SPEDWty( BIS OKs) (2.23a) 
and from Eq. 2.21c 
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Using the relation of Eq. 2.23a to substitute for [F] and tal invkge 


2.23b we get 
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* 
which proves that matrix [K]" of Eq. 2.23b is identical to [K] of Eq. 
iks3e; 
To demonstrate that the right hand side of Eq. 2.22 is equiv- 


alent to the right hand side of Eq. 1.3d and 1.4, call the inter-boundary 


* 
force, partition of Eq: 2.22 {Ry} » then 


zine a *4 2.23d 
{Ri} 7 {Ry} [F] ed {R. } ( 3 ) 


now substitute for [F] from Eq. 2.23a into Eq. 2.23d 
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* 
Vector {R,} of Eq. 2.23e is identical to {R,} ” of Eq. 1.3d. 


b 
For obvious practical considerations, submatrices eae fea 
[D], FKdr. and {Ry } ‘ will be obtained by modifying the coefficients of 
bas rand) s{Rierethesstep [D]T VEL ]i2 URAL Sime) 2.22 wilt be: re= 
stricted toi |i {R.} , since this is all that is needed until back- 
substitution begins. 
The relations required for backsubstitution can be obtained by 


carrying out the multiplications involved in Eq. 2.22. Thus, from the 


hirse row Of \E@. 2.22 


EOile Me tha) S40 H bl needa teed gen (Um (2. 24a) 
where 


* = 
{R, } ey {R. } (2.24b) 
Assuming that vector {U,} will be obtained from a solution of the 
inter-boundary system, vector {U, } , which represents the displacements 


of the internal nodes, can be obtained from Eq. 2.24a as 


ED tRy = Eel Ue (2.25) 
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Ceri (v.} (2.26) 


If the problem is not substructured, [K..], {U,} , and {R.} 
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default to [K], {U} , and {R} respectively, and Eqs. 2.20, and 2.26 
default to Eqs. 2.11 and 2.14a 


2.9 Columnwise Algorithms for Substructure Displacement Analysis 


In this section the algorithms necessary to carry out the 
partial release procedure described in Section 2.8 are derived. A 
number of algorithms can be derived for Eqs. 2.21, 2.25, and 2.26. Ifa 
skyline is imposed on the stiffness matrix, the resulting algorithms 
form the basis of a highly efficient equation solving package, suitable 
for substructure analysis. The derivation of the various algorithms, 
along with their final form after imposing the skyline on the stiffness 


matrix, follows. 


* 
2.9.1 Algorithm for [L]', [D], and {R,} 

Consider a 6x6 substructure stiffness matrix, where degrees of 
freedom 5, and 6 form the inter-boundary partition. Eq. 2.2la can be 
put in the form 
2d K Kits iK 
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Bay tae ‘a3 ! daa Kyg | (2-27) 


Expanding the fourth column will yield the set of equations 
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Koq = 944 bgy bay + do0 bao (2. 28b) 
K3q = 941 bay b31 * doo bao Lge + 433 by3 (2.28) 
K3q = 911 bay b31 + 4o0 Lago bao + 433 bg3 Lag + dag (2. 28d) 


Assuming that ae is known for i = 1,4, and j = 1,2, Eq. 2.28c can be 


solved for Ln3 as follows 


2 
Patypeeasenas,” “34° <2) Gea Las tie 232) 
Ln3 = G,/d34 (2.29b) 
and Eq. 2.28d can be solved for dag 
3 
day = Kaa ce Gog Las (2.29c) 


Generalizing Eqs. 2.29 yields the algorithm to determine foe and [D]. 
* 
Eqs. 2.24b can be put in the form {R. } = [L] {R. } , which 


is expanded for the above substructure as follows 
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P: Number of internal d.o.f. 
n: Total number of d.o.f. 
mm: Maximum of m;, m; 


Figure 2.7 Combined Algorithm for [L], [D], and {r;*} 
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Expanding the fourth row of Eq. 2.30 will yield 


(223) 


(2232) 


which when generalized, yields an algorithm for {R, } ; 
Notice that when Eqs. 2.29 and 2.32 are implemented on the 
computer Las will overwrite Koq> and dag will overwrite Kaa: These 


equations can be combined into the algorithm shown in Fig. 2.7. 


2.9.2 Algorithm for [F]! 
With Fist and [D] determined, as in Sect. 2.9.1, Eq. 2.21b can 


be expanded as follows, 
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where the kK, represent the elements of ba This notation is used 
because the elements of fri. overwrite the [KJ partition of the stiff- 


ness matrix. The expansion of the second column will yield 


= (2.34a) 
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Figure 2.8 Algorithm for [F] ™ 
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Kg = Rig 417/01 * Ko6 doo (2. 34b) 
K36 = Ki6 417 £31 + Kae doo b32 + Keg 433 (2. 34c) 
Kae > 6 411 bar * Ko6 doo Lago + K36 433 Lag + Kag daa (2. 34d) 
Eq. 2.34d can be solved for Kae if Kee. i = 1,3 are known. 

s. cw 

Kg = (Kyg - 2 Keg Keg Keg )/Kag (2.35) 
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Eq. 2.35b, when generalized, can be expressed by the algorithm shown in 
FAQs. .2.0 
* * 
2.9.3 Algorithm for [K] , and {Rt 
a re i 
With [F]” known and represented by [K], the components of [K] 


may be determined by expanding Eq. 2.2]1c 


K 55 K 56 Kes Keg 
* * 
N65 * 66 Kes Kee 
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- 111 Kyg> doo Kags 433 K36> Gna Kap] | Ko5 Kae 
K35 K36 
eae (2020) 
Component Ke can be expressed as 
* sr —_ — — — — — pe 
K cg = Kee - [day Kis Kig * tao Kos Kog + 433 Ky5 Kag + qq Kgs Kag] 
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Figure 2.9 The Combined Algorithm for [K]*, and {Re} 
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* 
With {R. } determined as in Section 2.9.1, and Ge deter- 


° ° . * 
mined in Section 2.9.2, {Ri} can be determined by expanding Eq. 2.23d 


as follows 
Rr) = Rts K 7 
Be | — | ys Kos Kg5 Kgs R | 
* a Paes et =o * 
R'6 Re K16 Kae Ko6 Kag R 5 
* 
R 3 
* 
R 4 (2.38) 
Component Re can be expressed as 
R™ * ok * a * +k * 
g ~ Re - Eye Ry + Kog Roa t Kgg Rog t Rag R a] (2.39a) 
or 
er 4 : me 
Rig = Re - 2 ae (2.39b) 


A combined algorithm representing Eqs. 2.37b, and 2.39b is 
* 
shown in Fig. 2.9. The inter-boundary stiffness matrix [K] of Eq. 2.22 
* 
and the associated force vector {R.} of Eq. 2.23e, have therefore been 


determined. 


2.9.4 Algorithms for Backsubstitution 


* . . 
Submatrices [kK] and {Rh} are the inter-boundary partitions. 
The assembly of these partitions from a number of substructures forms a 


higher level substructure unit in what has been called a hierarchy of 
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Figure 2.10 Algorithms for Backsubstitution 
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substructures. When the summit of the heirachy, referred to later as 
the ‘master system', has been formed, backsubstitution to obtain the 
final inter-boundary nodal displacements, and subsequently the internal 
nodal displacements, begins. 

Eqs. 2.25, and 2.26 have been expanded to the algorithm shown 
in Fig. 2.10, where the upper part of the algorithm evaluates VJ 
of Eq. 2.25, and the lower part of the algorithm is the backsubstitution 
OtieeGg., 2.26. 


2.9.5 The Equation Solving Package 


The algorithms developed above have been coded into a set of 
subroutines for solving substructured problems. The subroutines, along 
with their listing are described in Appendix C, and their use within the 


programs is described in Chapter 4. 


2.10 Efficiency of the Skyline Method 


The efficiency of this method relative to other methods is not 
readily determined in terms of the number of numerical operations since 
it is dependent ae the skyline of the sparse matrix, and on the partitions 
adopted. However, several aspects of the method indicate that a sub- 
stantial reduction of CPU time can be achieved. The saving in storage 
due to exclusion of zero wedges is obvious. The elimination of these 
wedges further eliminates the need to check for zero elements in the 
outer loops of the coding. Moreover, the method, as coded, automatically 
adjusts for the smaller vector when evaluating the scalar product of 


two vectors, whereas in banded algorithms, the vector that appears in the 
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outer loop governs the length of vector products for the next m operations. 
The fact that the elements required in evaluating the .vector 

products are accessed in sequential storage locations indicates that the 

technique will operate efficiently in virtual memory machines, and the 


Singly subscripted arrays may also be advantageous in some computers. 
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CHAPTER 3 - ASSEMBLY AND COORDINATE TRANSFORMATION 


3.1 Introduction to Assembly and Coordinate Transformation 


As shown in Fig. 3.1, a local reference frame is associated 
with each substructure. The position of a substructure in a higher 
level assembled system is defined by the inter-boundary nodal connec- 
tivity, and by the orientation of the local substructure reference frame 
relative to the higher level reference frame. The assembly of substruc- 
ture inter-boundary partitions into a higher level system must take care 
of both connectivity and coordinate transformation. 

Two assumptions were made to simplify the assembly process. 
The first, already mentioned in Chapter 2, is that the inter-boundary 
nodes of a substructure are grouped at the end of its local nodal num- 
bers. The second, is that there will be a fixed number of degrees of 
freedom per node, and that the local ordering of these degrees of free- 
dom will be constant throughout the problem. The second assumption 
implies that a stiffness matrix can be partitioned into a set of qxq 
Submatrices, where q is the number of degrees of freedom per node. Let 
these submatrices form the elements of a matrix [Q] equivalent to the 
stiffness matrix. By definition, component LQ] om is the qxq set of 
forces at node 2 due to a set of unit displacements applied separately 


at node m, when all other sets of displacements are zero. 


3.2 Coordinate Transformation 


Let vectors {U3 and {R)} of Eq. 2.18 be the displacement 


and force interboundary partitions of a substructure in its local 


reference frame, and let vectors {UJ and {Rp be the corresponding 
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Figure 3.1 Reference Frames 


Figure 3.2. Transformation Matrices 
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vectors in a higher level system reference frame. The orthogonal trans- 


formation between the two systems takes the form 
We Seana Eig (3.1a) 
* 
ORF Soler R: (3.1b) 


where [T,] and [T] are each composed of n? submatrices of order (qxq) 
and n is the number of inter-boundary nodes. These submatrices are null 
everywhere except those on the main diagonals of LT, ] and [T,], denoted 
respectively as [T,] and [T,] as illustrated in Fig. 3.2. 

Let u, v, and r be the local displacement components at a node 
in a plane frame problem, where u is the displacement component in the 
local X-direction, v is the displacement component in the local Y- 
direction, and r is an anti-clockwise rotation; and Up» Vi» and ry, be 
the corresponding components at the same node in a higher level system 
reference frame. Let Re Rie and Re be the force components at a node 
in the local reference frame, corresponding to degrees of freedom u, v, 
and r respectively, and let re ie and Ee be the force components at 
the same node in a higher level system.reference frame. The orthogonal 


transformations between the two sets of free vectors takes the form 
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where C is the cosine of angle 6, and S is the sine of 6, as illustrated 


in Fig. 3.3. Consequently, from Eqs. 3.1 and 3.2, 
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One 09 @ 


and the assembled transformation matrices for the inter-boundary nodes 


become 
Mig (it bee LT (3.4) 


The second row of Eq. 2.18 can be written in the form 
oS * 
OO Ry ee ie (3.5) 


Premultiplying Eq. 3.5 by ban we get 


AME IRC ht ee Ge ne! (3.6) 


where {F} is the vector of effective forces in the higher level ref- 
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erence frame. Since [K ] is symmetric, and [T] is orthogonal, the 
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Figure 3.3 Orientation of Reference Frames 
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Figure 3.4 Transformed Submatrix [Q],,, (A plane frame problem) 
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product (Ty! [kK] [T] is also symmetric, and represents the orthogonal 
transformation of [K'] from the local reference frame to the higher level 
reference Frame. Call the transformed stiffness matrix [Q]. If ia 
is a qxq nodal partition of Reds the corresponding partition of [Q], 


[Q],,,2 can be written as 
foie anny UK 11 (3.7) 


* 
Notice that [K lige and [Q] on are not symmetric if 2 #m. The explicit 


form of this transformed matrix LQ) om is shown in Fig. 3.4. 


3.3 Assembly 
Following the definition of [Q], the assembly of an inter- 


boundary partition in a higher level system matrix can be carried out 
in the same manner as the usual direct stiffness assembly. However, 
rather than assembling individual stiffness coefficients, qxq submatrices 
are assembled into qxq spaces in the higher level system stiffness 
matrix. In other words, instead of identifying the degrees of freedom 
of a substructure inter-boundary partition with the corresponding degrees 
of freedom of a higher level system, we need only identify the nodes of 
the inter-boundary partition with the nodes of the higher level system. 
This nodal connectivity information is very similar to the element nodal 
connectivity of simple direct stiffness techniques, and can be supplied 
either with every substructure, or with the higher level system data. 

The assembly scheme is shown in Fig. 3.5. Some difficulties 
arise from the fact that we are dealing only with upper triangles of the 


matrices. For example, submatrix [Q]o.4 for a substructure system, defined 
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by Eq. 3.7, may assemble into QR in the higher level system, where 
nodes 2 and m of the substructure identify with nodes k and j of the 
higher level system. Since we have only the upper triangle of the 
substructure stiffness matrix, m 2 eed sg 4 k, submatrix [Q],.; occupies 
a place inside the upper triangle of the higher level system matrix. If 
j < k, submatrix le is outside the upper triangle. But, since the 
higher level system matrix must be symmetric, then [Q), = 01, is 
assembled rather than [Q) 5: This procedure is shown in dotted lines in 
biG) (3...55 

Following the first assumption of Section 3.1, when the first 
inter-boundary node in a substructure is identified, the higher numbered 
nodes are all inter-boundary nodes, and can be taken in order. The 
assembly then can proceed across the columns of the inter-boundary 
partition of ren q columns at a time. 

In conclusion, it must be pointed out that the operations of 
transformation, according to Fig. 3.4, and assembly have been carried 
out simultaneously in order to reduce the operations needed to access 


the components of the qxq submatrices. 


3.4 The Skyline of a Higher Level Unit 


The assembly process described above controls the skyline of 
the higher level unit and access of [5] depends on the location of 
[55] which can be found from MAXA for the higher level unit. Estab- 


lishing this skyline can take several forms. This will be discussed 


briefly in the following. 
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3.4.1 The Direct Approach 


If the inter-boundary partition of a lower level unit-is 
assumed fully populated, it can be considered as an element, and the 
Process of computing and updating the column heights for the higher 
level unit would be the same as that described in Section 2.7.1. How- 
ever, this approach is simplistic, since the assumption of fully pop- 
ulated inter-boundary partitions will result in more storage space than 
is actually required, as well as a definite increase in the number of 
numerical operations performed during assembly, decomposition, and 
backsubstitution. This approach is the one primarily used in the pro- 
grams developed in Chapter 4, and forms the basis of the storage and CPU 
time analysis in Chapter 5. However, two alternative schemes to adjust 
the storage and access for the higher level unit to the ‘exact' skyline 


requirement have been developed, as indicated below. 


3.4.2 Skyline Modification 


Using the direct approach described above, it is expected that 
for some columns, the actual first non-zero component will be somewhere 
down the column from the location indicated by the direct higher level 
assembly procedure. Thus, if this can be checked after assembly, the 
array MAXA modified accordingly, and the components of the stiffness 
storage vector shifted forward as may be necessary (Fig. 3.6), the true 
skylines of the lower level units would have been taken into consider- 
ation in an indirect manner. This approach will result in a possible 
lower skyline for the higher level unit, and a saving in the number of 
numerical operations required for decomposition and backsubstitution, 


but not in storage requirements, since this modification must necessarily 
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take place after the assembly of the higher level unit stiffness matrix. 
This modification has been coded in subroutine MODAX1, see Appendix C, 
and the resulting saving in CPU time as applied to program MUSAPF, can 


be found in Table 5.2 denoted by ‘modification 1'. 


3.4.3 Predicting the Skyline 


In this Section a rigorous approach to establishing the skyline 
of a higher level unit, taking into consideration the skylines of the 
lower level units, is developed. The scheme is very similar to the 
assembly scheme described in Section 3.3 and Fig. 3.5. Utilizing the 
two assumptions mentioned in Section 3.1, the following algorithm may 
be established. 

For all the inter-boundary nodes (m) of the lower level unit 
taken in order, the following operations are performed. 

iP. For all submatrices [Q, ad of the lower level inter-boundary 
partition, starting at the main diagonal and ending at the 
level of partition or the skyline, whichever is lower, ident- 
ifty nodes k and j in the higher level system corresponding to 
nodes 2 and m respectively in the lower level system. 

2. ek E j, the proposed column: heights for the degrees of free- 
dom of node j are functions of the difference (j-k). If these 
column heights are greater than the current column heights, 
update the latter. 

3: If k > j, the proposed column heights for the degrees of free- 
dom of node k are functions of the difference (k-j). If these 
column heights are greater than the current column heights, 


update the latter. 
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This approach establishes a skyline which is somewhere between 
the skylines determined by the approaches of Sections 3.4.1 and ‘3.4.2. 
It has been coded in subroutine SKYPRD, see Appendix C. The resulting 
storage requirements and CPU time can be found in Tables 5.1, and 5.2 as 


‘modification 2'. 


3.5 External Boundary Conditions 


As stated in Section 3.1, there will be a fixed number of 
degrees of freedom per node, and the local ordering of these degrees of 
freedom will be constant through the problem. The method of imposing 
the external boundary conditions on the stiffness matrix must not violate 
this assumption. The method used in programs SISAPF and MUSAPF is to 
attach a spring with a very high stiffness relative to the structural 
stiffness coefficients, to the node where a degree of freedom is restrained 
in the direction of this degree of freedom. This method, in addition to 
its simplicity allows for specifying non-zero displacements. 

The coding of programs SISAPF and MUSAPF allows the user to 
impose the external boundary conditions at the basic substructure level 


and/or at a higher level of assembled substructures. 
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CHAPTER 4 - DEVELOPMENT OF PROGRAMS 


4.1 Introduction to Program Development 


Two substructure displacement analysis programs have been 
developed. Program SISAPF: Single-level Substructure Analysis for 
Plane Frames is a one level substructure analysis program, while 
MUSAPF: Multi-level Substructure Analysis for Plane Frames is a multi- 
level substructure analysis program. Both programs handle only plane 
frame type problems with no limitations on the size, number of substruc- 
tures, and number of cases of loading. 

The two programs use the equation solving package developed in 
Chapter 2 and described in Appendix C, as well as the coordinate trans- 
formation and assembly schemes described in Chapter 3. All subroutines 
developed are common to both programs, with the exception of the main 
executive routines and the input routines. A dry run facility is incor- 
porated in both programs to check the input data, and to determine the 
size of storage required for the stiffness matrices. 

The programs use a type of dynamic array allocation, proposed 
by McCormick [8]. This simple data manager identifies an array by its 
FORTRAN name, its length, and a logical number which designates the common 
block to which the array belongs. Further, the data manager can free 
the core storage assigned to the last defined arrays in any common block, 
or alternatively, it can free an entire common block. 

Data transfer to and from backing storage is done using unfor- 
matted read and write statements on sequential files. To achieve a high 
degree of flexibility in program MUSAPF, two IBM system routines, NOTE 


and POINT, are used to move the read and write pointers freely over the 
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sequential files. 


4.2 SISAPF: A Program for Single-level Substructure Analysis 


of Plane Frames 
This program performs a one level substructure displacement 
analysis of plane frames. It is written for large problems which possess 
a degree of regularity in structural properties. This Section describes 


the organization and flow of SISAPF. 


4.2.1 Concepts and General Description of SISAPF 


The program is composed of two main phases. The first phase 
is user controlled and involves the definition and formulation of the 
problem. The second phase is the solution and output phase. 

The definition of the problem has several aspects. A structure 
can be divided into a number of substructures, for which the inter- 
boundary connectivity and orientations can be considered global para- 
meters, and do, in fact, define the master system. Yet, some of these 
substructures may be identical in structural properties relative to a 
local reference frame. If these identical recurrent substructures have 
the same number of inter-boundary nodes, they may be completely defined 
by specifying the structural nionanbie of a representative structural 
unit (defining a substructure 'class'), together with the global sub- 
structure identification numbers of the individual substructures rep- 
resented by this class. 

The inter-boundary connectivity of these individual substructures 


describes the master system. This information is entered as array NPG. 
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Each column of this array holds the node numbers in the master system to 
which the substructure's interboundary nodes attach, starting with the 
first interboundary node of the substructure at the top of the column 
and taking the rest of these nodes in order. Table 4.1 describes array 
NPG for the master system of the example in Fig. 4.1. 

Within a class of substructures the same loads may be imposed 
On a number of substructures. Instead of entering the loads of each 
individual substructure separately, only the independent cases of load- 
ing within a class of substructures are entered in one array, together 
with a load case identification array, KSUB, which assigns to each sub- 
structure a number of local cases of loading equal to the number of 
global cases of loading. These loads are referenced to the local ref- 
erence frame. For example, the building shown in Fig. 4.1a can be 
divided into 5 substructures as shown in Fig. 4.1b. Substructures with 
global numbers 2, 3, and 4 have the same structural properties, and the 
same number of interboundary nodes. Thus they form a class of substruc- 
tures, say class number 2. Substructures number 1 and 5 form substruc- 
ture classes 1 and 3 respectively. This example has two cases of loading, 
shown in Fig. 4.1c, and it can be seen that substructure units 3 and 4 
have identical loads for global load case 1, while substructure unit 2 
has a different pattern. For global load case 2 all three substructures 
have identical loads. Thus inside substructure class number 2, there 
are three independent cases of loading that together describe all the 
global loads on all units inside this class. These three loading patterns 
can be entered in one array, together with the load case identification 


array KSUB shown in Table 4.2. 


Since the technique of defining the connectivity and loading 
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Figure 4.1b The Substructures 
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Figure 4.1c Global Load Cases on The Structure 


Figure 4.1 Details of Examnle for SISAPF 
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Array NPG 


Table 4.1 Array NPG for Master System 
of Example in Figure 4.1 


Local Substructure Number Be ay eo eee 
Global Substructure Number (NGSUB) ps toga ame 


Local Load Case Number for GLC. 


Local Load Case Number for GLC. 2 


Table 4.2 Arrays NGSUB and KSUB for Substructure 
Class 2 of Example in Figure 4.1 
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of a problem has now been described it is possible to return to a descrip- 
tion of the program. The MAIN program segment is composed of: a number 
of common statements; a number of statements where the lengths of the 
common blocks are specified and passed to the data manager; a call 
statement to the main executive routine MAINMG; and, an END card. In 
this way the lengths of the common blocks can be exactly calculated 
prior to run time, and passed to the program through this small program 
segment which is compiled at run time. The program can be used to 
predict the required lengths of the common blocks, by employing the dry 
run facility described in Appendix A. 

The executive routine MAINMG calls all other subroutines, 
which read the data, formulate and solve the problem, and output the 


final results. 


4.2.2 Logic Flow of SISAPF 


The sequence of events coded into the main executive routine 
MAINMG of program SISAPF is as follows. Details of the required input 


are given in Appendix A. 


1. Read the problem control parameters (Sub: INPUT1). 

on If the problem is unsubstructured, go to 9. 

oF Read the master system connectivity information, as well as sub- 
structure orientations (Sub: INPUT2). 

4. Read the master system external boundary conditions, if any (Sub: 
BOUND) . 

2) Form the array of column heights for the master system stiffness 


matrix, by looping over the connectivity data supplied in step 3 


(Sub: COLHT):: 
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18. 


19%: 


20. 


5G: 


Form the addressing array for the diagonal components of the master 
stiffness matrix in the master system stiffness vector (Sub: 
ADDRES). ) 

If this is a dry run, go to 9. 

Reserve and clear space for master system stiffness and load vectors. 
Set ICLSUB = 1, and start to loop over substructure classes, in 
order to read the data, form, reduce, and assemble the stiffness 
and load vectors. 

Read substructure class ICLSUB control parameters (Sub: INPUT3). 
Read nodal geometry, member material and cross section properties, 
member connectivity data, and arrays KSUB, and NGSUB for substruc- 
ture class ICLSUB (Sub: INPUT4). 

Read external boundary conditions, if any (Sub: BOUND). 

Read nodal loads, and/or prescribed displacements (Sub: JLOAD). 
Read member loads, and, if this is not a dry run, compute the 
corresponding nodal loads (Sub: MLOADS). 

Form the column height array for the substructure class ICLSUB 
stiffness matrix (Sub: COLHT). 

Form the diagonal component addressing array for the substructure 
class ICLSUB stiffness vector (Sub: ADDRES). 

Ltetnis. 1s a dry run; go to 23: 

Form and assemble member stiffnesses into the substructure class 
ICLSUB stiffness vector (Sub: STIFF). 

Add external boundary conditions if any to the substructure class 
ICLSUB stiffness vector (Sub: BOUND2). 

If required, decompose the stiffness partition [K..] and reduce the 


load partition {R, } (Sub: EQSBST) according to the algorithm of 
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Fig. 2.7. If the problem is unsubstructured, go to 28. 

Reduce the stiffness partition [K.] to tae and form the inter- 
boundary partitions [kK ] and Rr according to the A tuortenn 
in Figs. 2.8 and 2.9 (Sub: EQFT and Sub: EQKBB). 

For each substructure that belongs to class ICLSUB, assemble 


* 
partitions [K ] and {R into the master stiffness and load 


s 
vectors, with the appropriate coordinate transformations (Sub: 
ASSEMB) . 

If the total number of classes of substructures equals 1, go to 31. 
Pr Chiss1s-aedry run, "gO “to 26: 

Store array pointers, control parameters, and all relevant arrays 
Oresubstructure "class TCESUB on Pile =|: 

Free common blocks SUBIA, SUBRA, and SUBSR. 

If ICLSUB equals the total number of classes of substructures, go 

to 31. Otherwise, go to 10 to start on the next class of substruc- 
tures. 

This is a branch to an unsubstructured problem, reached from step 
20. Backsubstitute from the decomposed stiffness matrix into the 
reduced load vector, to obtain the nodal displacements (Sub: BKSB1). 
Output displacements and member end forces:( Sub: sDISPI® andy sub: 
STRESS)” 

Return to MAIN. End of unsubstructured problem. 

1f-this ts.a dry runs “go tov4Z. 

Add external boundary conditions, if any, to master system stiffness 
vector (Sub: BOUND2). 


Solve for master system nodal displacements (Sub: EQSBST and Sub: 


BKSB1). 
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34. Start the backsubstitution and output control loop. Rewind FILE 1, 
and set ICLSUB = 1. 

35. If the total number of classes of substructures equals 1, go LOe3/* 

36. Read array pointers, control parameters, and all relevant arrays of 
substructure class ICLSUB from FILE 1. 

37. Start to loop over all substructures that belong to class ICLSUB. 
The global identification numbers can be obtained from the array 
NGSUB. For each substructure, steps 38 to 40 will be performed. 

38. Copy the inter-boundary nodal displacements from the master system 
solution vector into a dummy array BA, referenced to the local 
reference frame of this substructure. Also copy partition ik 
from the substructure class loading array into BA. The local load 
case number corresponding to any global load case can be obtained 
from array KSUB (Sub: RESUB). 

39. Backsubstitute, if required, from the decomposed substructure class 
stiffness matrix into the dummy vector BA, to obtain the substructure 
internal nodal displacements (Sub: BKSB2, and Sub: BKSB1). 

40. For each global load case, output the nodal displacements, and member 
end forces for this substructure. 

41. Set IGLSUB = ICLSUB + 1... If ICLSUB is less than or equal, to the 
total number of classes of substructures, go to 36 to start on the 
next class of substructres. 


42. Return to MAIN. End of problem. 


The function of subroutines INPUT], INPUT2, INPUT3, and 


INPUT4 is described in steps 1, 3, 10 and 11 respectively. Subroutines 
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COLHT, ADDRES, EQSBST, EQKBB, EQFT, BKSB1, and BKSB2 form the equation 
solving package, and are listed in Appendix C. Subroutines BOUND, JLOAD, 
MLOADS, STIFF, BOUND2, ASSEMB, DISPL, STRESS, and RESUB are commen to 
programs SISAPF and MUSAPF, and are listed in Appendix C. The data 
transfer to and from FILE 1 is carried out by four short subroutines, 
contained in Appendix C. Finally the data manager is also common to 


both programs and is listed in Appendix C. 


4.3 MUSAPF: A Program for Multi-level Substructure Analysis 


of Plane Frames 
This program is a more advanced version of program SISAPF, that 
deals with large complex problems that can be decomposed into a heirarchy 
of substructures. In its present form, it can handle only plane frame 
type problems. It has no limitations on the number of substructures or 
the number of levels of the heirarchy, and also has no limitations on the 
number of cases of loading. This Section describes the organization and 


flow of MUSAPF. 


4.3.1 Conceptual Development 


The concept of a heirarchy of substructures is outlined " 
Chapter 1. However, an illustrative avennie is useful to describe how 
the concept has been implemented. The illustrative structure shown in 
Fig. 4.2a can be partitioned into 9 basic structural units (Fig. 4.2b). 
Two overlapping numbering systems will be used to designate the various 
substructures. The first, a ‘global numbering' system identifies each 


substructure as a physical component. In Fig. 4.2b, global units 1 and 
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9 have the same structural properties, and hence they form a class of 
substructures for which the structural configuration will be designated 
by ‘independent basic unit' number 1. This independent basic unite will 
be denoted as IBU (1). Global unit numbers 2, 4, 6, and 8 have identical 
structural properties, which will be designated by IBU (2). In the same 
‘manner global unit numbers 3, 5, and 7 can be represented by IBU (3). 
Global units 1 and 2 can be assembled along their inter-boundary nodes, 
as shown in Fig. 4.2c. This new assembled unit is a higher level sub- 
structure unit. Give this unit global number 10, and since it is an 
independent structural entity, designate it IAU (4), which stands for 
‘independent assembled unit' number 4. In the same manner global unit 
pairs 3 and 4, 5 and 6, 7 and 8, can be assembled to form global units 
11, 12, and 13, respectively. These three new higher level substructure 
units have identical structural properties, having been assembled from 
the inter-boundary systems of 3 identical pairs of independent units. 
All three may be represented by IAU (5). Further, global units 10, 11, 
12, 13, and 9 may be assembled together along the interboundary nodes 
(Fig. 4.2d), giving a still higher level substructure unit, which has 
global number 14, and may be designated by IAU (6). This last assem- 
blage is unique, and forms the top of the hierarchy of substructures, 
call it the master system. 

The loads can be prescribed physically at the independent 
basic unit level. For example, for the loading situation shown in Fig. 
4.3a. IBU (2) has three cases of loading as shown in Fig. 4.3c. These 
cases of loading are independent loading patterns, which are pooled in 


one array, and which can represent all loading combinations on global 
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units 2, 4, 6, and 8. Load combinations affecting independent assembled 
units can be represented by inter-boundary forces arising from the 
possible load combinations on lower level units. The loading eee at 
all levels can be traced without difficulty from Fig. 4.3. 

The control arrays for an IAU are INDSUB, IRCSUB, IBNSUB, NPG, 
LCSUB, and ORINT. These arrays are described in detail in Appendix B. 
However, for IAU (5) of the example, Table 4.3 shows the first three of 
the arrays, Table 4.4 shows array NPG (the connectivity array), and 
Table 4.5 shows array LCSUB. 

Once the master system has been assembled, its nodal displace- 
ments can be obtained, and the backsubstitution process is begun to 
solve successively for the nodal displacements of all substructures. 
This process can follow several paths, but the main feature is that a 
unique solution vector must be obtained for each global unit of the 
heirarchy. The problem is to minimize the data transfer to, and from 
backing storage. Since it is necessary for the user to specify the way 
the substructures are to be combined into higher level global units, it 
has been decided to let the backsubstitution process also be controlled 
by the user, since after going through the former process the user has 


all the information required to specify the latter process. Output can 


be initiated once a basic unit has been reached. Table 4.6 illustrates 


a user prepared complete backsubstitution process for the above example. 
Row 2 of this table holds the global unit number of a higher level unit 


for which a solution vector exists, while row 3 has the global number of 


a lower level unit for which a solution vector is sought. Row 4 supplies 


the independent unit number representing the lower level unit. Row 5 


gives the position of lower level unit in the control arrays of the 
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Figure 4.2b Basic Substructure Units 
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Figure 4.2c Assembled Substructure Units 


Figure 4.2d The Master System 


Figure 4.2 A Hierarchy of Substructures 
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Figure 4.3a The Global Loading Pattern 


Figure 4.3b Loads on IBU(1) 
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Figure 4.3c Loads on IBU(2) 


Figure 4.3d Loads on IBU(3) 
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Figure 4.3e Loading Pattern on IAU(4) 
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Figure 4.3f Loading Pattern on [AU(5) 


ale ae ee ee 


\ IOUROERRCRRRER! ISURRURORERR0R! (QURCRRRCR0008) 


= (ERURECRCRERERI 


Constituent Unit Numbers: 4 5 5 5 1 
Load Case Numbers for L.C.1: 1 1 2 2 1 


Figure 4.3g Loading Pattern on IAU(6) 


Figure 4.3 Description of Hierarchial Loading 
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Table 4.5 Array LCSUB for IAU(5) 


INDSUB 
IRCSUB 
IBNSUB 
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independent unit representing the higher level unit. For example, IAU 
(5) which represents unit 13 occupies position 4 in the control arrays 
of IAU (6) which represents unit 14, the 'master system'. This last 
parameter is necessary for the identification of the proper load cases 
of the independent lower level unit forming the global loading scheme of 
the lower level unit. The execution of Table 4.6 is carried out column 
by column. For example, executing the first column, the inter-boundary 
displacements of unit 10 are copied from the solution vector of unit 14, 
into a dummy array BA, referenced to the local reference system of unit 
10. The proper load cases of IAU (4), representing the loads on unit 10 
are obtained using the parameter of the last row, and the corresponding 
internal loading partitions are copied from the local load array of IAU 
(4) into the dummy array BA. Next, the stiffness matrix of IAU (4) is 
retrieved from backing storage, and backsubstitution to compute the 
internal nodal displacements of global unit 10 is carried out. As 
global unit 10 is not a basic substructure unit, array BA, the solution 


vector, will be stored on a file, and execution of the second column 


starts. 


4.3.2 General Description of MUSAPF 


The program uses the same data managing and equation solving 
packages as program SISAPF, as well as the same data transfer, formula- 
tion, and output subroutines. The MAIN program segment is similar to 
that of SISAPF. The difference between the two programs is that MUSAPF 
uses two types of substructures (basic units and assembled units), and 
also backsubstitution is a user controlled process in MUSAPF, as described 


in Section 4.3.1. The higher level control arrays are ORINT for orienta- 
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tion, and NPG for connectivity of the constituent units, as well as 
LCSUB which identifies the local load cases on a constituent unit which 
defines its effect on the higher level unit. The input description is 


found in Appendix B, and the listing in Appendix C. 


4.3.3 Logic Flow of MUSAPF 


The sequence of events coded into the main executive routine 


MAINMG of program MUSAPF is as follows. 


1. Read problem control parameters (Sub: INPUT1). 

2. Read backsubstitution control arrays (Sub: INPUT2). 

Sia Set ICLSUB = 1, and start to loop over the independent basic sub- 
structures units (IBU's), to read data, and form, reduce, and store 
the IBU stiffness and load vectors. 

4. Read IBU(ICLSUB) control parameters (Sub: INPUT3). 

5. Read IBU(ICLSUB) nodal geometry, member material and cross section 
properties, and member connectivity data (Sub: INPUT 4). 

6. Read IBU(ICLSUB) external boundary conditions, if any (Sub: BOUND). 

7. Read IBU(ICLSUB) nodal loads, and/or prescribed displacements, if 
any (Sub: JLOAD). 

8. Read IBU(ICLSUB) member loads, and if this is not a dry run, compute 
member end forces (Sub: MLOADS). 

9. Form the column height array for this IBU stiffness matrix (Sub: 
COLHT). 

10. Form the diagonal component addressing array for this IBU stiffness 
vector (Sub: ADDRES). 


Wa olf this 1s aedryi run, go to 18. 
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We: 


19. 


20. 


ae 


Reserve and clear space for this IBU stiffness vector. 

Form and assemble member stiffnesses into the IBU stiffness vector 
(Sub: STIFF). | 

Add external boundary conditions, if any, to this vector (Sub: 
BOUND2). 

If required, decompose the stiffness partition [K;.] and reduce 
load partition {R. } according to algorithm of Fig. 2.7 (Sub: 
EQSBST). 

If this is an unsubstructured problem, go to 19. Otherwise, reduce 
the stiffness partition [K..] to the form cry! according to the 
algorithm of Fig. 2.8, and form the inter-boundary partitions eee 
and mR according to the algorithm of Fig. 2.9 (Sub: EQFT, and 
Sub: EQKBB). 

Store the present values of the 'write' pointers of FILE 1 and 2 in 
arrays IRTRVI and IRTRV2, respectively. Then store the array 
pointers, control parameters, and all relevant arrays of IBU(ICLSUB) 
OneELEES land 2, 

Set ICLSUB = ICLSUB + 1, and if it is greater than the total number 
of independent basic units, go to 22. If not, go to step 4 to 
start on the next IBU. , 

This is a branch to an unsubstructured problem. Backsubstitute 
from the decomposed stiffness matrix, into the reduced load vector, 
to obtain the nodal displacements (Sub: BKSB1). 

Output the nodal displacements and member end forces (Sub: DISPL, 
and Sub: STRESS). 


Return to MAIN. End of unsubstructured problem. 
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Start to loop over independent higher level units. These are 
assemblies of more than one IBU or IAU unit. Their numbering 
system is a continuation of the numbering of the IBU's. 

Read IAU(ICLSUB) control parameters (Sub: INPUTS). 

Read the constituent units' identification numbers, connectivity 
and load case identification arrays (Sub: INPUT6). These constituent 
units may be IBU's or IAU's provided that they have been previously 
defined. 

Read IAU(ICLSUB) external boundary conditions, if any (Sub: BOUND). 
Form the column height array for the IAU(ICLSUB) stiffness matrix, 
by looping over the connectivity data supplied in step 23 (Sub: 
COLHT). 

Form the diagonal component addressing array for the IAU(ICLSUB) 
stiffness vector (Sub: ADDRES). 

Ri thisitis.fa, any runs? go vto (36. 

Reserve and clear space for the IAU(ICLSUB) stiffness and load 
vectors. 

Assemble the constituent units', inter-boundary stiffness and load 
partitions, into the IAU(ICLSUB) stiffness and load vectors with 
the appropriate transformations (Sub: ASSEMB). The constituent 
unit matrices can be obtained from FILE 2, for which the ‘read’ 
pointers are obtained from array IRTRV2. 

Add external boundary conditions, if any, to the IAU(ICLSUB) stiff- 
ness vector (Sub: BOUND2). 

If required, decompose the stiffness partition [K..] and reduce the 


load partition {R.} according to the algorithm of Fig. 2.7 (Sub: 


EQSBST). 
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If this is the master system, go to 38. 

Reduce the stiffness partition [K. J to the form au and form the 
interboundary partitions [Rage and (R} according to the algo- 
rithms of Figs. 2.8 and 2.9, respectively (Sub: EQFT, and Sub: 
EQKBB). 

Store the present values of the 'write' pointers for FILES 1 and 2 
in arrays IRTRV1, and IRTRV2 respectively, and store the array 
pointers, control parameters, and all relevant arrays of IAU(ICLSUB) 
ine RILES. 1. and: 2. 

Free common blocks from this IAU's arrays. 

If this is a dry run, and if this is the master system, go to 41. 
Set ICLSUB = ICLSUB + 1. Go to 23 to read the data, and defint the 
next IAU. 

Obtain the solution vector for the master system (Sub: BKSB1), and 
if required, store this vector in FILE 3. 

Start the backsubstitution and output control loop. This process 
is an execution of a table similar to Table 4.6, which is supplied 
in step 2. If required, solution vectors can be stored on FILE 3, 
for which the write pointers are stored in array IRTRV3. 


End of problem. Return to MAIN. 


The input subroutines INPUT] to INPUT6 are briefly described 


in steps 1, 2, 4, 5, 23, and 24, and they are not to be confused with 


Subroutines INPUT] to INPUT4 of program SISAPF. 
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CHAPTER 5 - APPLICATION 


5.1 Description of Structure 
A flat slab high rise structure [7] is analysed, using SISAPF 


and MUSAPF, both as a one unit structure and as a partitioned structure 
in several ways. An analysis of the execution times for several schemes 
used to partition the building is attempted. 

In Figs. 5.la and 5.1b, the equivalent frame [6] and a half 
plan of the building are shown. The columns and shear wall are of 
constant section and the slab is seven and a half inches thick. All 
sections are assumed uncracked. The link members are of high axial 
stiffness to simulate rigid floor diaphragms. The shear wall is approxi- 
mated by a vertical member of the correct stiffness, and horizontal 
rigid projecting beams are used to simulate wide column to beam con- 
nections [14]. 

The structure is analysed under a lateral wind load of 20 psf 
intensity in the Y direction, see Fig. 5.1b. These loads are applied to 
the joints along the left side of the lumped frames and on the left side 


of the shear wall-frame. 


5.2 Partitioning Schemes 


Five different schemes have been used to partition the struc- 
ture. The first scheme is an analysis of the structure as one unit and 
will be referred to as Example 1 when run on MUSAPF, and Example 6 when 
run on SISAPF. The second scheme is a single level substructure scheme, 
Fig. 5.2, and will be referred to as Example 2 when run on MUSAPF, and 


as Example 7 when run on SISAPF. In this scheme, there are 15 basic 
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Figure 5.1 


Figure 5.1a Elevation of frame 


Figure 5.1b Typical Floor Plan of Structure 


Details of Problem Structure 
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Figure 5.2 Substructure Scheme for Examples 2 and 7 
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Scheme for Example 3 


Figure 5.3 Substructure 
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Figure 5.5 Substructure Scheme for Example 5 
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(b) Member Numbering 
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Figure 5.6 Details of IBU(1) for Examples 2 and 7 
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Figure 5.7 Details of IBU(3) for Examples 2 and 7 
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Figure 5.8 Details of IBU(2) for Examples 2 and 7 
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Figure 5.9 Details of 1AU(4) or the Master System for Examples 2 and 7 
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substructure units divided into three classes, and one assembled unit 
which is the master system. The details of these units are shown in 
Figs. 5.6, 5.7 and 5.8, the master system is shown in Fig. 5.9, and the 
input files for both programs are in Appendix D. The third, fourth, and 
fifth schemes are multi-level substructure schemes, see Figs. 5.3, 5.4, 
and 5.5, respectively, for which the details have not been provided. 
These have been run on MUSAPF. 

The analysis of all examples yielded identical displacements 
and member end forces for all nodes and members up to 6 decimal places 
which proves the accuracy of the equation solver, and the fact that the 


multi-level substructure scheme used, does not affect the accuracy. 


5.3 Size and Core Storage 


One of the main advantages of the substructure method is the 
saving in core space required to process a problem. Table 5.1 gives the 
control parameters, the master system stiffness size, and the maximum 
size of core space for each example in normal length words. It can be 
seen from Examples 1, 2, 3, and 4, that as the number of levels of sub- 
structures increases, the core space required decreases. The maximum 
core space required for any substructure scheme is between 35% and 47% 
of the size required for the structure as one unit. As indicated in 
Section 3.4, this advantage could be further enhanced, if the skylines 
of the lower level matrices were considered when computing the skylines 
of the higher level unit stiffness matrices. The discussion above is 
based on the ‘direct assembly' approach, described in Section 3.4.1. 

The effects of modifying the skylines, using the assembly schemes described 


in Sections 3.4.2 and 3.4.3 are discussed in Section 5.6. 
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5.4 CPU Time 

Programs MUSAPF and SISAPF use the IBM system subroutine TIME 
to output the CPU time at every stage of execution. The values obtained 
depend to a certain degree on random machine conditions. However, a 
comparison is possible, if all examples are run in one session, taking 
care to specify the correct file sizes. 

The stages of execution are, mainly: the reading and formulation, 
the decomposition, the storage on peripheral devices, and the backsub- 
stitution and output. The CPU time taken by each stage is shown in Table 
5.2 for all examples, as well as, the total CPU time. 

It can be seen that all substructured examples, with the 
exception of Example 4, take less CPU time than the unsubstructured 
examples. The CPU time saving should increase if the element type 
requires large formulation time. The formulation time is an average 
of 15% of the total time for the substructured examples as opposed to 
50% for the decomposition stage. The time required for backsubstitution 
is an average of 35% of the total time for the substructured problems. 
The backsubstitution time increases as the number of levels of substruc- 
tures increases, whereas the decomposition time depends mainly on the 


size and number of the stiffness matrices. 


5.5 Substructuring, and Nodal Numbering 


In all examples, the nodal numbering, for both the basic level 
substructures and the assembled units, has been chosen carefully. The 
governing factor was the element connectivity for the former, and the 


constituent unit connectivity for the latter. Over a number of experi- 
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mental runs not included in those examples, it has been observed that 
numbering the master system vettically instead of horizontally as in 
Fig. 5.9, raises the decomposition CPU time by as much as 500%. 
Partitioning a structure is better carried out parallel to the 
Shorter side. This way the master system will have fewer nodes, and any 
particular substructure will span fewer interboundary nodes, resulting 
in higher efficiency. The reason is the fact that the master system 
nodes are necessarily inter-boundary nodes in some or all substructure 
units at all levels. The repeated assembly of, and backsubstitution 
through these nodes consumes more time as their number increases. 
Numbering the inter-boundary nodes of any substructure, or the 
nodes of the master system parallel to the shorter side lowers the 
skyline of this particular system in the regions [K; J. and [K J: 


Lowering the skyline means less storage and less CPU time. 


5.6 Higher Level Unit Modified Skylines 


The modifications to skylines of higher level units discussed 
in Section 3.4, have been applied experimentally to program MUSAPF. The 
approach characterized by checking the columns of an assembled higher 
level unit stiffness matrix for the actual first non-zero component, is 
referred to as ‘modification 1', while the approach characterized by 
predicting the skyline of a higher level unit taking into consideration 
the skylines of the lower level units, is referred to as ‘modification 
ok 

Table 5.1 shows the change in stiffness storage vector length 
for modifications 1 and 2, for Examples 2 to 5. It must be noted that 


while modification 1 indicates a smaller storage requirement, this is 
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not true in practice, since the modification takes place after the 
assembly. However the storage saving for modification 2 is real. 

Table 5.2 compares the total CPU time consumed for Examples 2 
to 5, for both modifications with the unmodified scheme. For 'modifi- 
cation 1', the saving varies between 14.4% for Example 4 to 3% for 
Example 5, while for ‘modification 2', the saving varies between -0.4% 
for Example 5 to 12.3% for Example 4. Two trends can be observed; the 
first being the increased saving in CPU time as the number of levels of 
substructures increases, while the second is that as the nodal numbering 
approaches an optimum state, no significant change is observed. 

‘Modification 2' is more elaborate, and has a rigorous theoret- 
ical background compared to ‘modification 1'. However, the latter 
appears to be more economical for several reasons. The number of numer- 
ical operations involved in 'modification 1' can never be greater than 
that of ‘modification 2', and while the latter requires a certain amount 
of data retrieval, the former does not. Also ‘modification 1' gives a 
better upper bound for the skylines as can be seen from Table 5.1, since 
it checks for the actual first non-zero component in a column, while 
‘modification 2' anticipates this component. The above reasoning is 
supported by the fact that 'modification 1' gives consistently less CPU 
time, and should be recommended as a permanent addition to programs 


MUSAPF, and SISAPF. 
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CHAPTER 6 - SUMMARY AND CONCLUSIONS 


Two plane frame substructure analysis programs, SISAPF and 
MUSAPF have been developed. SISAPF is based on a single-level substruc- 
ture scheme, and MUSAPF is based on a multi-level substructure scheme. 
For this purpose, an equation solving package based on the skyline tech- 
nique utilizing the concept of 'Cholesky' decomposition, as well as an 
assembly and a coordinate transformation scheme have been developed. 

From studies on a sample structure, it has been observed that 
using a substructure scheme the core space requirements can be reduced 
to 35% of that required for the same structure as one unit. A saving in 
CPU time amounting to 15% can be easily achieved. The saving in CPU 
time should be greater if more complicated structural elements are used 
Since the saving in formulation time averages 37%. While the saving in 
decomposition time averages 26%, a rise in backsubstitution CPU time is 
observed and increases as the number of, and number of levels of, sub- 
structures increases. 

Partitioning and inter-boundary node numbering are found to 
give best results when both are carried out parallel to the shorter side 
of the structure or substructure. 

As they stand, both programs accept only one type of element. 
However, both can be developed further to accommodate structures where 
the number of nodes per element, and the number of degrees of freedom 
per node can be varied. Generalization to accommodate such changes 
would require modification of the assembly and coordinate transformation 


schemes. 


The basic assembly scheme assumes that the inter-boundary 
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partitions of lower level substructure units are fully populated when 
assembling a higher level unit. This scheme has been improved by recog- 
nizing the skylines of the lower level unit stiffness matrices when 
assembling the higher level units, and developing two modified assembly 


schemes which take this into account. 


code Patel eer LET SRB ean 


\ salle ,. ct 
“gover Vd OSV Maga ea sah 9280 RS 


pat esr! ay eT ae re! ¢! ver 


Weed bet 7 oe Yt eT OVE LH 
diy 
ae 
H 
Sag 
ae 
<0 
7s § 
i 
= 
, 
i 
a 
4 
_ 
i 
7 
= > 


10. 


Te 


86. 


REFERENCES 


Bathe, K.J., and Wilson, E.L., 'NONSAP - A Nonlinear Structural 
Analysis Program', Nuclear Engineering and Design 29 (1974), 
pp. 266-293. 


Bathe, K.J., and Wilson, E.L., "Numerical Methods in Finite Element 
Analysis', Prentice Hall Inc., Englewood Cliffs, New Jersey, 
1976. 


Felippa, C.A., ‘Solution of Linear Equations with Skyline-Stored 
Symmetric Matrix', Lockheed Palo Alto Research Laboratory, 
Palo Alto, California, 1973. 


Furuike, T., ‘Computerized Multi-level Substructuring Analysis', 
Computers and Structures, Vol. 2, pp. 1063-1073, Pergamon 
Press, Great Britain, 1973. 


Hrabock, M.M., 'The Method of Substructures and its Application to 
Structural Analysis', M.Sc. Thesis, University of Saskatchewan, 
Regina, Saskatchewan, 1975. 


Khan, R.F., and Sbarounis, J.A., ‘Interaction of Shear Walls and 
Frames', Journal of the Structural Division, ASCE, Vol. 90, 
No. ST3, June, 1964. 


Mathews, G.S., 'The Second-Order Analysis and Design of Reinforced 
Concrete Frames', M.Sc. Thesis, University of Alberta, Edmonton, 
Alberta, 1977. 


McCormick, J.M. 'Data Handling in FORTRAN Programs', Committee on 
Computers and Information Systems, American Society of Civil 
Engineers, National Structural Engineering Convention, 

975. 


Meyer, C., ‘Solution of Linear Equations - State of the Art’, 
Journal of the Structural Division, ASCE, Vol. 101, No. ST4, 
April 1975, pp. 869-890. 


Murray, D.W., 'PLFRAM: Program for Plane Frame Direct Stiffness 
Analysis', University of Alberta, Edmonton, Alberta, 1975. 


Przemieniecki, J.S., ‘Theory of Matrix Structural Analysis', McGraw- 
Hill, New York, New York, 1968. 


Rosen, R., and Rubenstein, M.F., ‘Substructure Analysis by Matrix 
Decomposition’, Journal of the Structural Division, ASCE, 
Vol. 963. No. S13. March 719705 pp.) 663-670. 


Schrem, E., ‘Computer Implementation of the Finite Element Procedure’, 
Numerical and Computer Methods in Structural Mechanics, 
Fenves, S.J., Perrone, N., Robinson, A.R., and Schnobrich, W.C., 
Ed., Academic Press, Inc., New York, New York, 1973, pp. 79- 
hel 


ier wt te 70a bt Lito A. sep angi! a 
| he lval one gn entand a 


‘ae 

if Rodel xo ined 

. wo eT TO bows haa 2on) eae 

fs bd mu fy r ? wha te! 45 
ba nn ri4 a fa 
ave 

Oi meh cet il 9 2gue fc a ~ ia Gti 
mA OTERO TD i So. FO Bh iy 


j honk Pf. Krein smh 9 APY Ledbic 4 
Lest ae Tri vt? 13 WI Mb ss yt end; 


ait STON Vane Foo tee nl 


i? fo TORR, lore vad PA: vide 


NONI AD pol gall “bok 4 i oe c.9) tat 
NOFA aL ie ees aa aa} . ieee ae i 

sms seaTs. * ii 

| é' fwrickh’ » “unto? nsersion ne: 


‘ivid “30 -¢aSteao2 We saver <4 a hie eae 9 ; 
NIRV? LelhyseroInsy: PRT: de . wt .cmanntoa?, | 
: i} Ae ve : p 


290 chy Sy OPED .~ coms Yuertil ys G fart 3 Tere * eed) 9, “anys! 
Ee OM VM Toe, STR (oar eid Ve he Sut 4 Lane iv" 
( y > “G4 Goes. f INGA, av, 


akan 3282 Foal enter” gitel® vot nats mantis aby a 


OVO! Bagh on oserong . wyedth WH NSN a . Af eghs 


na 
werrdot 2 bay onl. '24gesgaherxt tal! vf y “sgt! Pama f St ieatngtnas Rites 
: a | GR Aa Pe fa aan 2) 


st eae, RNa | oak Sud outs ‘Me soot enPaitysi ai? Meas afta 208 


» 2024 pNETE IVT are ad 70 ah | es oad 
: c zO\8 LAG, i0F 4 , Wee aa a 
‘sy 


hae 


. an, Snamats win t4 we hy pine 4 ‘vamTane utd! > ee Ef a. : 

a 9 29nBi sel Ceui pede ih abuntiny Pasay hee Tag taped re oe 

2 bs) Tongs: GAs sh aeERtaGM. | Me ‘oe ebb’ » COW, iF 
) a a ota peeet aie wait =F avo ¥. 30h « CATHODE IA ys 

“oe 


ven ans } 
7 se | ] 7 ' 4 ® + 
Da on : —— : — 7 ay 


Bie 


14. Schwaighofer, J., and Microys, H.F., ‘Analysis of Shear Walls Using 
Standard Computer Program’, ACI Journal, December 1969. 


15. Williams, F.W., ‘Comparison Between Sparse Stiffness Matrix and 
Substructure Methods', International Journal for Numerical 
Methods in Engineering, Vol. 5, 1973, pp. 383-394. 
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APPENDIX A_- USER'S MANUAL FOR SISAPF 


A.1 Master System Data Cards 


A.1.1 Heading Card (20A4) 


One card which contains any title for the problem. 


A.1.2 Master System Control Card (8145 F 12.0, 14) 
4 So PAZ 16 20 24 28 32 44 48 
NSUB : Total number of substructures. 
NCLSUB : Total number of classes of substructures. 
MXND : Maximum number of inter-boundary nodes in any substructure. 
NUMNST : Total number of nodes in the master system. 
NGLC : Total number of global cases of loading. 
IDOF : Number of degrees of freedom per node = 3. 
NODE : Number of nodes per element = 2. 
MBEL : Total number of external boundary elements attached to the 


master system, if any. 

UNIT : Conversion factor for length units, e.g. 12.0 if section and 
material properties are entered in 'inch' units, while the 
nodal geometry is entered in 'foot' units. 

IDRY self, (sents isa, dryerun: 


If 0, this is a production run. 
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A.1.3 Substructure Orientation Cards (8F10.0) 
If NSUB = 1 on card A.1.2, this group is to be omitted. This 
group of cards holds the orientation of the substructure local frames of 
reference relative to the master system reference frame. Enter as many 


cards as necessary to describe all substructures (8 substructure 


units/card). 
10 20 30 80 
ORINT(I) : Local reference frame orientation relative to master system 


reference frame, measured in degrees (positive anti-clockwise). 


A.1.4 Master System Connectivity Data Cards (2014) 


If NSUB = 1 on card A.1.2, this group is omitted. The group 
describes the connectivity of each substructure with the master system, 
and forms a table of width NSUB, and length (MXND + 1), where each 
column of the table holds the information for one substructure. Sub- 
structures must be ordered according to the global substructure identi- 


fication number. An example of this table is shown in Table 4.1. 


--- NPL(J) 


--- NPG(1,J) 


J : The global number of the substructure. 


NPL (J) - Number of inter-boundary nodes in this substructure. 
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NPG(K,J) : Node number in the master system to which node K of sub- 
structure J attaches, counting from the first inter- 


boundary node in this substructure. 


NOTE: If there are more than 20 substructures, this table must be 
divided into a number of tables, entered successively, with the restric- 


tion that each table must contain (MXND + 1) cards. 


A.1.5 Master System External Boundary Element Cards . (3 Tari 220)) 


One card per external boundary element. If MBEL = 0 on card 
A.1.2, this group is to be omitted. A boundary element is defined as a 
| spring with very high stiffness attached to a node. Each ‘element' can 
restrain one or more of the degrees of freedom at a node. Several 
elements with different stiffnesses in different directions can be 


attached to one node. 


E ran wone) [AT 


: Number of the boundary element. 

MNPB(N) : Number of the master system node to which this element 
attaches. 

MKODE(N) : A three digit number of the form ijk corresponding to degrees 
of ee u, Vv, and r (where r is the nodal rotation), 
respectively. If any of these is 1, the corresponding 
degree(s) of freedom will be restrained; if any of these is 
0, the corresponding degree(s) of freedom is (are) unrestrained. 

AEB(N) : The spring stiffness of this element, if other than 107°. 


The units should be consistent with the material properties. 
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A.2 Substructure Class Data Cards 
This group will be repeated for every class of substructures. 
The types of data cards required for each class of substructure consists 
of the following: 
2.1 Substructure Class Control Card 
2.2 Nodal Geometry Cards 


2.3 Member Data Cards 


A. 
A. 
A. 
A.2.4 Substructure Load Case Identification Cards 
A.2.5 External Boundary Element Cards 

A.2.6 Joint Load/Displacement Cards 

A.2.7 Member Loads 


These cards are described in detail below. 


A.2.1 Substructure Class Control Card (814) 
Ae Orc. -uihO 20 24 28 32 
IC : Number of this class of substructure. 
NJ : Total number of nodes for this class. 
NE : Total number of elements (i.e. members) for this class. 
NLC : Total number of independent cases of loading including zero 


cases of loading, for this class. 


NSUBCL : Total number of substructures that belong to this class. 
NEBEL : Total number of external boundary elements for this class. 
NFIBN : Number of the first inter-boundary node for this class. 
NIBNS : Total number of inter-boundary nodes for this class. 


NOTE: NJ = NFIBN + NIBNS - 1 
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A.2.2 Nodal Geometry Cards CR 2 oO ee) 


One card/node, unless automatic generation of nodal data is 


used. 


AN X(N) | ¥(N) 


N : node number. 

X(N) : X-coordinate of the node. 

Y(N) : Y-coordinate of the node. 

INC : If non zero, automatic nodal generation will be initiated. 


The generated nodes will have numbers (NOLD + K*INC), where 
NOLD is the nodal number on the preceding card, and K is 
a positive integer. Generation terminates when (NOLD + 
K*INC) = N. The coordinates of the generated nodes are 
linearly interpolated between node NOLD and node N._ INC 


should be positive. 


NOTE: The first nodal geometry card must have INC = 0, and the last 


must have N = Nd. 


A.2.3 Member Data Cards 


A.2.3.1 Member Property Default Card (SE 12.0) 
KZ 24 36 
ADEF : Default value for member area. 
RDEF - Default value for moment of inertia of member. 


YDEF - Default value for Young's modulus of member material. 
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Notice that if UNIT = 12.0 on card A.1.2, ADEF, RDEF, and YDEF should be 
entered in in.*, in.*, k/in.* respectively. If UNIT = 1.0 any consistent 


set of units may be used. 


A.2.3.2 Member Property and Connectivity Cards (7A Sed2.0) 


One card per member, unless automatic member data generation 


is used. 


4 8 2 16 20 24 28 40 52 64 
| NOD(T.M) NOD(J.M) MKODE(M) | INC | INCI INCd | AREA(M) | RI(M) | yMop(M) 


M : Member number. (NOTE: The last member card must have M = NE) 

NOD(I,M) : End I node number. 

NOD(J,M) : End J node number. 

MKODE(M) : 0, if member is continuous at both ends. 

Le AE member is hinged at end I only. 
2, if member is hinged at end J only. 
3, if member is hinged at both ends. 

INC : If non zero, automatic generation of member data will be 
initiated. The new members will have numbers (MOLD + K*INC), 
where MOLD is the member number on the preceding card, and 
K is a positive integer. These members will have the cross 
section and material properties of the member on the card 
initiating the generation (i.e. the present card). The end 
nodes of the new members will be (NOD(I,MOLD) + K*INCI) and 
(NOD(J,MOLD) + K*INCJ) for nodes at end I, J respectively. 
INC must be zero for the first member card of the substructure, 
and zero or positive otherwise. 

INCI : Incrementation value for nodal number of end I. 


INCJ : Incrementation value for nodal number of end J. 
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AREA(M) : Area of member section, if other than ADEF. 
RI(M) : Moment of inertia of member section, if other than RDEF. 


YMOD(M) =: Young's modulus of member material, if other than YDEF. 


NOTE: Units of AREA, RI, YMOD are the same as for ADEF, RDEF, and YDEF. 


A.2.3.3 Echo Check Flag Card (14) 
A 
K : If 1, the complete nodal geometry and member data will be 


contained in the output. 
If 0, the output will not contain an echo check of the 


completed data. 


A.2.4 Substructure Load Case Identification Cards (2014) 

This group identifies the global numbers of the substructures 
which belong to this class and the local cases of loading on each which 
form the global loading pattern for all global load cases. The group is 
composed of a table the first row of which holds the global identifica- 
tion numbers of the substructures starting in the second tabular column. 


The rest of the rows (cards) each hold the local load case numbers which 


define a global load case. 


4 8 12 
le] NGSUB(1) |NGSUB(2) 
) 


K |KSUB(K, 1) | KSUB(K, 2) 


16 76 80 
--- (NGSUB(J), J=1, NSUBCL) --- 


--- ((KSUB(K,J), J=1, NSUBCL), K=1, NGLC) --- 
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J : Subscript to span the total number of individual substruc- 
tures*oip this: class. 

NGSUB(J) : Global identification number of this substructure. 

K : Global load case number. 

KSUB(K,J) : The local load case number that when imposed on global sub- 
structure number NGSUB(J), forms a part of global load case 
number K. These local substructure load cases will be 


specified in Sections. A.2.6 and A.2.7. 


NOTE: If NSUBCL is greater than 19, this table must be divided into 


several tables, each composed of (NGLC + 1) cards. 


A.2.5 External Boundary Elements Cards (Si aawercrO) 


If NEBEL = 0 on card A.2.1, this group is to be omitted. The 
group is composed of NEBEL cards, one card per boundary element. See 


Section A.1.5 for a definition of the external boundary elements. 


4 8 We 24 
nwo foe [ees 
N : Number of the external boundary element. 
NPB(N) : Number of the substructure node to which this element is 
attached. 
KODE(N) : A three digit number of the form ijk, corresponding to 


degrees of freedom u, v, and r, respectively. If any of 
these is 1, the corresponding degree of freedom will be 
restrained, otherwise the digit must be 0. 


BES(N) : The spring stiffness of this element, if other than 10°°. 
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A.2.6 Joint Load/Displacement Cards 


This group consists of a number of cards per loaded node, 
unless automatic load generation is used. The first card in a node 
group is a nodal identification card. The rest of the cards for this 


node each describe a local non-zero case of loading. 


A.2.6.1 Nodal Identification Card (714) 
4 8 12 16 20 24. — 78 
[Jee fnese [owe [owe wets [ne 
N : Node number (these need not be in order). 
NMLC : Number of non-zero cases of loading at this node. 
IDESP : If non-zero, displacements will be specified for one or more 


cases of loading. 

(NBEL(I), : The boundary element numbers attached to this node, in case 

I=1,3) non-zero displacements are specified. These three numbers 
correspond to degrees of freedon u, v, and r, respectively. 
In case this node is restrained in 3 directions with one 
element, and non-zero displacements are specified only in 
2 directions, then this boundary element number appears 
only twice in the corresponding locations. 

INC : If non-zero, automatic loading generation will be initiated. 
Nodes with numbers (NOLD + K*INC), where NOLD is the node of 
the preceding nodal group, and K is a positive integer, wil] 
be assigned the same loads as On the node initiating the 
generation (N) for all local cases of loading. N must be 


greater than NOLD for successful load generation. 
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A.2.6.2 Nodal Load Cards (214s USF W203) 
One card/non-zero case of loading. Number of cards of this 


group is NMLC, as input on card A.2.6.1. 


IgG : Number of local load case (need not be in order). 

JKODE(LC) : A three digit number of the form ijk, the digits corresponding 
respectively to u, v, and r. If any digit is 1, the corres- 
ponding u, v or r will be interpreted as a displacement in 
the direction of the corresponding NBEL of card A.2.6.1, 
otherwise it must be zero. 

U(N,NLC) : Load in the direction of the X-axis. 

V(N,NLC) : Load in the direction of the Y-axis. 


R(N,NLC) : Moment in the X-Y plane, positive anti-clockwise. 


NOTE: If UNIT = 12, and YMOD is in K/in.*, u, and v must be in kips and 
r in K.ft. Also notice that the u, v, and r arrays do not exist in the 
program. The corresponding values are entered directly in the local 


load array, B. 


A.2.6.3 Termination Card 


One blank card at the end of the ‘joint load/displacement card 


group'. 


A.2.7 Member Load Specification Cards 


This group consists of a number of cards for each element on 


which non zero loads are applied, unless automatic load generation is 
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used. Any number of loads can be applied on an element in one or more 
load cases. Load generation is carried out for one load for one case of 
loading at a time, and it can be of two types. The first type generates 
loads on elements which do not have the same end conditions, nor the 
same orientation or length. The second type involves members which have 
the same orientation and length, but not necessarily the same end condi- 
tions. In cases other than when automatic load generation is used, the 
cards may be entered in any order, since any one card defines completely 
the member number, the local load case number, the load, and its type, 


position, and orientation with respect to the member. 


A.2.7.1 Member Load Cards (AL4. O2EV2702"4 E1020) 
One card per member per load per load case, unless automatic 
load generation is used. 
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M : Member number. 
EC : Local load case number to which this load belongs. 
INC : If non zero, automatic load generation will be initiated for 


load case LC. The load on the card initiating the generation, 
or its end effects, will be assigned to members with numbers 
(MOLD + L*DABS(INC)), where MOLD is the member number on the 
preceding card, which must have the same LC value, and L is 

a positive integer. 

If INC is positive, the generation is of the first type, and 
if INC is negative, the generation will be of the second type 


described in Section A.2./7. 
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CL 
CN 


QI 
AI 
QJ 
AJ 


NOTE: 


: If 1, the load is concentrated 


If 2, the load is of uniform intensity 


If 3, the load is trapezoidal. 


: Projection along the member of a vector in direction of load. 


: Projection normal to the member of a vector in direction of 


load. 


: Magnitude of load at distance (AI*length) from end I. 
: Fraction of length from end I to beginning of load. 
: Magnitude of load at distance (AJ*length) from end I. 


: Fraction of length from end I to end of load. 


If the load is concentrated, QJ and AJ may be omitted. If the 


load is uniform, QJ may be omitted. Units of load are in K, or K/ft., 


1 UNE =< 122.0; 


A.2.7.2 Termination Card 


One blank card at the end of the group of member load cards. 


(The program returns to A.2.1 for input data for the next substructure). 


100. 


Loi ry 


ié T| in a 7 WAT ve ivew TAY i % eT: ' aa 


| Wes : pan 
fa wr boli mh evga wae ew 


au 
> Ss 


vie 


ter Wee rears “EA); os 
rint Fe ae rr oe ane ‘ot ’ 
i J ‘ 


ia Wo i ier! "OA nee 
oe A : 


008 4) ate Letedeiy 


std attri 94.4e OR tie OC ah. : 
he ee _ ey nih, us 


hada = 


; beat’ in 
22 aD e iy! *ecie 70) alah B aif 2 ie sii Lie aula a. 
laura Aaee Jen ang et we gqnt vr ad OF Srcatow ime RO. 


A.3 Common Block Description 


The size of the problem is passed to the data manager through 
the MAIN segment, as indicated in Section 4.2.1. The calculation of the 
sizes of the different common blocks should be carried out prior to run 
time and values assigned to array ICOM in segment MAIN according to the 
following formulas. The values of NWK and NWA may be set to zero for 
the dry run described in Section A.3.6. The dry run outputs the values 
of these variables to allow recomputation of common block sizes prior to 


the production run. 


A.3.1 Size of Common Block MASTRA 
ICOM(1) = NSUB + IDOF*NUMNST*NGLC + MBEL + NWK 


A.3.2 Size of Common Block MASTIA 
ICOM(2) = NAUB*(MXND + 1) + IDOF*(2*NUMNST + MXND) 


+ 2*(MBEL + NEBEL) + | 


A.3.3 Size of Common Block SUBRA 
ICOM(3) = IDOF*NLC*(NJ + NODE*NE) + 3*NE + 2*NJ + NWA 


A.3.4 Size of Common Block SUBIA 
ICOM(4) = NE*(NODE + 1) + NSUBCL*(NGLC + 1) + 2*IDOF*NJ 


+ IDOF*NODE + 1 


A.3.5 Size of Common Block SUBSR 
ICOM(5) = IDOF*(NODE*NE + NGLC*NJ) 
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A.3.6 The Dry Run Facility 


This option is implemented to check the input data, and to 
calculate the exact sizes of the stiffness storage vectors NWK and NWA. 
The option is exercised by specifying IDRY as indicated in Section 
A.1.2. During this run the sizes of the common blocks MASTRA, and SUBRA 
are input with zero values of NWK and NWA, where NWK is the size of the 
master system stiffness vector, and NWA is the size of a class of sub- 
structure stiffness vector. The dry run will give these values, which 
are then used to update the corresponding common block sizes prior to 
the production run. For this purpose the upper bound of ICOM(3), ICOM(4) 


and ICOM(5) values over all classes of substructure should be used. 
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APPENDIX B_  - USER'S MANUAL FOR MUSAPF 


B.1 Problem Control Cards 


B.1.1 Heading Card (20A4) 


One card which contains any title for the problem. 


B.1.2 Problem Control Parameter Card (Sia Elz. 0.714) 
4). 8 12 16 20 24 28 32 44 48 
IDOF | NODE | UNIT | IDRY 
NSUB : Total number of basic substructure units. 
NINDSB : Total number of independent basic units 
NTOTAS : Total number of substructures and substructure assemblies 


including the master system. 
NINDAS : Total number of independent units including the basic 


independent units, the assembled units, and the master 


system. 
NGLC : Total number of global cases of loading. 
NMBKSB : Number of backsubstitution steps. 
IDOF : Number of degrees of freedom per node = 3. 
NODE : Number of nodes per element = 2. 
UNIT : Conversion factor for length, e.g. 12.0 if the section and 


material properties are in 'inch' units, while nodal geometry 
1S) ine feeu.e uni ts. 
IDRY If ale thises a dry runs 


1f0;) Ehtse isa production run. 
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B.2 Backsubstitution Control Cards (2014) 

This group of cards describes the backsubstitution and is 
similar to Table 4.6. The table is composed of four rows (cards) and 
NMBKSB columns. If NMBKSB > 20 the table must be divided into several 
tables each consisting of four cards. If NTOTAS = 1, the problem is 


unsubstructured and this group must be omitted. 


4x] 
--- INFIBK(I) --- 
--- INF2BK(1) --- 
--- INF3BK(I) --- 


76 


I : Subscript denoting the number of the backsubstitution step. 


INFIBK(I) : Global number of a higher level unit for which a solution 
vector has been determined, e.g. INFIBK(1) = NTOTAS. 

INF2BK(I) : Global number of a lower level unit for which a solution 
vector is sought. It must be < NTOTAS. 

INF3BK(I) : Number of the independent unit which represents the lower 
level unit. It must be < NINDAS. 

INF4BK(I) : Local number of the lower level unit in the higher level 


unit control arrays. 


105: 


tera) Br 

#4 NitG) NOT ds ba Hi phe psy 
te ASG | eae WN) Tuy availa ah 
yor atat babivib esd: daums oer Ser Ans 


: es, rs : 
y a. 
andes] SATIN *i «ih A 


a 


on gag iba Ae ety ie faa ‘agi 
eee "tebe a) waaay 

sawot of exicagrpry Ho FN Wy ay -Hiaiining oid te. ‘at Paes ; 
me 

‘baked ot foun a stir Dowel. free 


oust eadbtd eat At aban! Reet Set a 19; Noein? heawd sql Tyeant’ | 
| 
| ont Horan0% ‘yin |} ie 
ae su ey 
rah * 
es ; at 
in di 
oe ; 
: 7 oa 
is at ? 
c 4, iF 


106. 


B.3 Independent Basic Unit Data Cards 


This group will be repeated for every independent basic unit. 


The types of data cards required for each unit are as follows: 


B.S. | 
bee ears 
ePRCH IS) 
B.3.4 
Bess 


B.3. 1 


Independent Basic Unit Control Card 
Nodal Geometry Cards 

Member Data Cards 

External Boundary Element Cards 


Loading Data Cards 


These cards are described in detail below. 


Independent Basic Unit Control Card (714) 


4 


Le 


NJ 
NE 
NLC 


NEBEL 
NFIBN 
NIBNS 


ool 


16 20 


24 28 
NEBEL | NFIBN | NIBNS 


: Independent basic unit number. The units must be entered 


in order. 


: Total number of nodes for this unit. 
: Total number of elements for this unit. 


: Total number of independent cases of loading for this unit. 


If a zero loading on the unit occurs in the global loading 
pattern, the zero loading must be considered a case of 


loading. Thus, if this unit is not loaded at all NLC = 1. 


: Total number of external boundary elements, if any. 
: Number of first inter-boundary node. 


: Total number of inter-boundary nodes for this unit. 


B.3.2 Nodal Geometry Cards (TAS 2612405 04) 


One card per node unless automatic nodal generation is used. 
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4 16 28 32 
n | X(N) y(n) Inc | 


The definition of these variables can be found in Section A.2.2. 
B.3.3 Member Data Cards 


B.3.3.1 Member Property Default Card VSrZa0)) 


12 24 36 


The definition of these symbols can be found in Section A.2.3.1. 


B.3.3.2 Member Property and Connectivity Cards C7 T4e 3riz. 0) 


One card per member unless automatic member data generation 
is used 


4 52 64 


The definition of these variables can be found in Section A.2.3.2 


B.3.3.3%, Echo.Check. Flag Card (14) 
4 
K : If 1, the complete nodal geometry and member data is output. 


If 0, the output will not contain an echo check of the 


completed data. 


B.3.4 External Boundary Element Cards (314. SR1Z30) 


One card per external boundary element. If NEBEL = 0 on card 
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8 12 16 20 24 28 40 
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B.3.1, this group is to be omitted. The definition of a boundary element 
can be found in Section A.1.5. 
4 8 12 24 
nN | NPB(N)| KODE(N) | BES(N) 
The definition of these variables can be found in Section 


ey aeeod, 


B.3.5 Load Data Cards 


B.3.5.1 Load Flag Card (214) 
4 8 
JLFLAG| MLFLAG 
JLFLAG : If zero, no joint loads will be prescribed on this unit. 


If non-zero, joint loads will be prescribed on this unit. 
MLFLAG : If zero, no member loads will be prescribed on this unit. 


If non-zero, member loads will be prescribed on this unit. 


B.3.5.2 Joint Load/Displacement Cards 


This group consists of a number of cards per loaded node, 
unless automatic load generation is used. The first card in a nodal 
load group is a node identification card. The rest of the cards ina 
nodal group, each describe a local non-zero case of loading. If 


JLFLAG = O on card B.3.5.1, this group is omitted. 


B.3.5.2.1 Node Identification Card (714) 


One card at the head of a nodal load group. 
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4 8 12 16 20 A a TAS: 
NMLC} IDESP | NBEL(1) | NBEL(2) | NBEL(3) | INC 
The definition of these variables can be found in Section 


Psce0el. 


B.3.5.2.2 Nodal Load Cards (214, 3F12.0) 
One card per non-zero case of loading on the node identified 
Oneecard’ Bb. Sa5e2. 1% 
4 8 20 32 44 
Lc | KODE (Lc) U(N,NLC) |V(N,NLC) | R(N,NLC) 
The definition of these variables can be found in Section 


WeesOnce 


B.3.5.2.3 Termination Card 


One blank at the end of the joint load/displacement card 


group. 


B.3.5.3 Member Load Specification 


This group consists of a number of cards for each member on 
which non-zero loads are applied, unless automatic load generation is 
used. Any number of loads can be applied on an element in one or more 
load cases. Load generation is carried out for one load for one local 
case of loading at a time, and it can be of two types. The first type 
generates loads on members which do not have the same end conditions, 
nor the same orientation or length. The second type involves members 
which have the same orientation and length, but not necessarily the same 


end conditions. In cases other than when automatic load generation is 
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used, the cards may be entered in any order, since any one card defines 
completely the member number, the local load case number, the load, and 
its type, position, and orientation with respect to the member. If 


MLFLAG = 0 on card B.3.5.1, this group is omitted. 


B.3.5.3.1 Member Load Cards (ATARI 250564 B00) 
One card per member per load per load case, unless automatic 
load generation is used. 
16 28 40 50 60 70 80 
nfucQne x To foufox area as 
The definition of these variables can be found in Section 


Ase</ G1 


B.3.5.3.2 Termination Card 

One blank card at the end of the group of member load cards. 
(The program returns to B.3.1 for input data for the next independent 
basic unit, unless IC on card B.3.1 for this group = NINDSB, in which 


case it proceeds to B.4 for the independent assembled unit data cards). 
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B.4 Independent Assembled Unit Data Cards 


unit. 


This group must be repeated for each independent assembled 


They must be in order, ending with the master system. If the 


problem is unsubstructured (NSUB = NTOTAS = 1), this group is to be 


omitted. 


The types of cards required for each independent assembled 


unit consists of the following: 


B.4.] 


Independent Assembled Unit Control Card 


B.4.2 Constituent Unit Information Cards 


B.4.3 External Boundary Element Cards. 


B. 


These cards are described in detail below. 


Independent Assembled Unit Control Card (914) 


eve NEBELH | NFIBNH ine at Oe feu ea | 


NJH 
NLCH 
NEBELH 
NFIBNH 
NIBNSH 
NCU 
MXND 


: Independent assembled unit number. These are a continuation 


of the IC numbers on card B.3.1, and are entered in order 
ending with the master system which must have ICH = NINDAS. 
If any of the constituent units is itself an assembled unit, 
it must have a lower ICH number, i.e. it must have already 


been read and defined. 


: Total number of nodes for this unit. 

: Total number of independent load cases for this unit. 

: Total number of external boundary elements for this unit. 
: Number of the first inter-boundary node for this unit. 

: Total number of inter-boundary nodes in this unit. 

: Total number of constituent units for this assembled unit. 


: Maximum number of inter-boundary nodes in any of the constituent 
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units for this assembled unit. 

IFLAG : This variable is only meaningful for the IAU which constitutes 
the master system. It should be non-zero only if the master 
system solution vector is to be stored for a backsubstitution 
that requires this vector and occurs after backsubstitution 


for lower level units has commenced. 


B.4.2 Constituent Unit Information Cards 


B.4.2.1 Connectivity and Load Case Identification Cards (2014) 
This group forms a table of width NCU and length (MXND + 3 + 


NLCH), where these variables are defined on card B.4.1, where every 
column holds the necessary connectivity information for one constituent 
unit. If the number of constituent units NCU is greater than 20, the 
table must be split into several tables with the restriction that each 
table must have (MXND + 3 + NLCH) cards. 
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INDSUB(I) : 


IRCSUB(I) : 


IBNSUB(I) : 


NPG(J,1) 


ECSUB(Ks.D):: 
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: A subscript which denotes the position or local number of a 


constituent unit in the assembled unit ICH. 

Identification number of the independent unit which represents 
the constituent unit I. INDSUB(I) < ICH. 

If zero, it means that the arrays of independent unit 
INDSUB(I) are already in core. This flag helps to eliminate 
the time needed to retrieve these arrays from backing storage, 
in case this independent unit represents more than one 
constituent unit. Such constituent units should occupy 
adjacent columns of this table with the first having 

IRCSUB(I) = 1. 

Total number of inter-boundary nodes for constituent unit 


INDSUB(I). 


: Number of the assembled unit ICH node to which node J of 


constituent unit INDSUB(I) attaches, counting from the first 


inter-boundary node in the lower level system. 


: A subscript which denotes the load case for the assembled 


unit ICH. 

The identification number of the independent load case out 
of system INDSUB(I), which when applied on constituent unit 
I forms a part of the loading scheme for load case K on 


assembled unit ICH. 


B.4.2.2 Constituent Unit Orientation Cards (8F10.0) 


Any number of cards to describe the orientation of the local 


frames of reference of the constituent units. (Eight units per card). 
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ORINT(1) --- ORINT(1) Solace oNCU 
I : A subscript which denotes the position or local number of a 
constituent unit in the assembled unit ICH control arrays. 
ORINT(I) : The orientation of the frame of reference of constituent 
unit I relative to the frame of reference of the assembled 


unit ICH, measured in degrees (positive counter-clockwise). 


B.4.3 External Boundary Element Cards (SI45¢F 1260) 


One card per external boundary element. If NEBELH = 0 on card 
B.4.1, this group is to be omitted. The definition of an external 


boundary element can be found in Section A.1.5. 
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The definition of these variables can be found in Section A.2.5. (The 
program now returns to B.4.1 for the input data of the next assembled 


unit, unless ICH = NINDAS, i.e. the last unit). 
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B.5 Common Block Description 


The size of the problem is passed to the data manager through 
the MAIN segment, as indicated in Section 4.3.1. The calculation of the 
sizes of the different common blocks should be carried out prior to run 
time, and values assigned to array ICOM in segment MAIN according to the 
following formulas. Where more than one formula is given, the largest ~ 
requirement governs. The values of NWK, and NWA may be set to zero for 
the dry run described in Section B.5.6. The dry run outputs the values 
of these variables to allow the recomputation of common block sizes 


prior to the production run. 


B.5.1 Size of Common Block PROBIA 
ICOM(1) = 4*NMBKSB + 3*NTOTAS ¢ 1 


B.5.2 Size of Common Block RA] 
ICOM(2) = 2*NJ + 3*NE + IDOF*NODE*NLC*NE + NEBEL 


B.5.3 Size of Common Block IAI] 
ICOM(3) = NE*(NODE + 1) + 2*NEBEL + NLC 
or ICOM(3) = NCU*(MXND + NLCH + 3) + 2*NEBELH 


B.5.4 Size of Common Block RAZ 
ICOM(4) = NCU + NJH*IDOF*NLCH + NWK + NJ*IDOF*NLC + NWA 
or ICOM(4) = NCU + NJH*IDOF*NGLC + NWA + NJ*IDOF*(NLC + NGLC) 


+ NE*NODE*IDOF. 


B.5.5 Size of Common Block IA2 
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ICOM(5) = IDOF*(NOD + 1 + 2*NJH) 
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or ICOM(5) = IDOF*(2*NJH + NJ + MXND) 
or ICOM(5) = IDOF*NJ + NE*(NODE + 1) 
or ICOM(5) = NCU*(MXND + NLC + NGLC) 


Beso. ou ihe Dry Russ Facility 

This option is implemented to check the input data, and to 
calculate the exact sizes of the stiffness storage vectors NWK and NWA. 
The option is exercized by specifying IDRY = 1 as indicated in Section 
B.1.2. During this run the size of the common block RA2 is calculated 
with zero values for NWK and NWA, where NWK is the size of the stiffness 
storage vector for an assembled unit, and NWA is the size of the stiff- 
ness storage vector for a lower level unit. The dry run will output 
these values, which are then used to update the common block RA2 size 
prior to the production run. For these values the upper bounds for al] 
the common block sizes over all possible configurations of the involved 


variables should be used. 
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APPENDIX C 
PROGRAM DESCRIPTIONS AND LISTINGS 
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APPENDIX C - PROGRAM DESCRIPTIONS AND LISTINGS 


C.1_ Program SISAPF 
Program SISAPF consists of the following parts. 


1. MAIN 
2. MAINMG, the main executive subroutine. 
3 The input package 
INPUT], INPUT2, INPUT3, and INPUT4 
4. The data managing package 
ISPAC, LOCOM, REMOV, REMOV2, and BLOCK DATA. 
5. The data storage and retrieval package 
CLEAR, ICLEAR, RTRV1, RTRV2, STORE], and STORE2. 
6. The equation solving package 
ADDRES; COLHT, -EQSBST; EQFT, EQKBB,, BKSB1,; and BKSB2. 
7. The formulation and output package 
ASSEMB, BOUND, BOUND2, DISPL, JLOAD, MLOADS, RESUB, 
STIFES and: STRESS: 
The listing for the first three packages follows. The listing 
for packages 4 to 7 can be found in Sections C.3, C.4, €.5, and C.6 


respectively. 
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C.2 Program MUSAPF 
Program MUSAPF consists of the following parts. 


1. MAIN 
2. MAINMG, the main executive subroutine. 
3. The input package 
INPUT], INPUT2, INPUT3, INPUT4, INPUT5, INPUT6. 
4. The data managing package 
ISPAC, LOCOM, REMOV, REMOV2, and BLOCK DATA. 
5. The data storage and retrieval package 
CLEAR, ICLEAR, RTRV1, RTRV2, STORE], STORE2. 
6. The equation solving package 
ADDRES, COLHT, EQSBST, EQFT, EQKBB, BKSB1, BKSB2. 
7. The formulation and output package 
ASSEMB, BOUND, BOUND2, DISPL, JLOAD, LOADID, MLLOADS, 
RESUB,) STIFF, and STRESS. 

The listing for the first three parts follows. The listing 
for parts «4,25, 6, and 7. can be found in;sections1€.3, 1C:4,°C.5,; and 
C.6, respectively. 

The modifications discussed in Sections 3.4.2, and 3.4.3 have 
been implemented by adding subroutine MODAX1 to the equation solving 
package for modification 1, and by subroutine SKYPRD to the same package 
for modification 2. The listing and description of both subroutines can 
be found in Section C.5. The corresponding changes to the main executive 


subroutine MAINMG are not shown herein. 
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C.3 The Data Managing Package 


This package is a modification of a simple manager proposed 
by McCormick [8]. It is composed of functions ISPAC and LOCOM, sub- 
routines REMOV, and REMOV2, and a BLOCK DATA segment which initializes 
the common block DIMCOM. 

The function ISPAC is called from subroutine MAINMG to define 
an array in one of five common blocks. This function then calls function 
LOCOM to check if this array is currently defined in the same common 
block. If it is not, then ISPAC enters the name of the function in 
the name directory (NAMES) corresponding to the required common block. 
ISPAC then adds the length of this array to the current length of this 
common block, and checks if the updated length of the common block 
exceeds the maximum specified in MAIN. Finally it returns with the 
position of the first element of this array in the common block. Sub- 
routine REMOV is called from MAINMG to remove any of the last defined 
arrays in any common block from the name directory (NAMES), and the 
pointer directory (IPT). Subroutine REMOV2 is called from MAINMG to 
free an entire common block by initializing the corresponding column 
in NAMES and IPT. The arguments of these segments are as follows 

ISPAC (NAME, LEN, K) 

LOCOM (NAME, K) 

REMOV (NAME, K) 

REMOV2 (kK) 
where, NAME is the name of an array to be defined, LEN is the array 


length, and K is a logical number to designate the required common block. 


The listing of this package follows. 
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C.4 Data Storage and Retrieval Package 

This package is composed of six subroutines. Subroutines 
CLEAR, and ICLEAR initialize real and integer arrays, respectively. 
Subroutines RTRV] and RTRV2 retrieve respectively integer and real 
arrays from backing storage, while subroutines STORE], and STORE2, store 
respectively integer and real array in backing storage. The arguments 
of these subroutines are as follows 

CLEAR (ARRAY, LEN) 

ICLEAR (IARRAY, LEN) 

RTRV1 (IN, IARRAY, LEN) 

RTRV2 (IN, ARRAY, LEN) 

STORE] (10, IARRAY, LEN) 

STORE2 (10, ARRAY, LEN) 
where, ARRAY is the first element of a real array to be initialized, 
retrieved, or stored, IARRAY is the first element of an integer array 
to be initialized, retrieved, or stored, LEN is the length of an array, 
IN,IO is a logical number which designates a sequential file. The listing 


of this package follows. 
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C.5 The Equation Solving Package 


C.5.1 Description of Package 


This package is based on the algorithms derived in Chapter 2 
for substructured columnwise decomposition and backsubstitution. It is 
composed of nine segments for which the description and listing follows. 


The arguments are described in Section C.5.3. 


C.5.2 Description of Subroutines 
(1) ADDRES (MAXA, MHT, NN, ID, NWA) 


This is a simplified version of subroutine ADRESS [2]. It 
computes the addresses of the diagonal components of a stiffness matrix 
stored columnwise in a vector, as well as the total length of this 
vector. 

(2) COLHT (NB, NODES, ND, ID, ME, MHT, NP, LM) 

This is a modified version of subroutine COLHT [2]. It forms 
and updates the active column heights of a stiffness matrix, see Section 
2.7.1. It is called per structural element which forms a part of a 
basic substructure unit, or per substructure unit which forms a part of 
a higher level assemblage. 

(3) EQSBST (A, B, MAXA, NP, ID, NLC) 

This subroutine is a modified version of subroutine COLSOL 
[2]. It performs a columnwise decomposition of the internal degrees of 
freedom of a substructure stiffness matrix and reduces the right hand 
side up to the same level. It is based on the algorithm shown in Fig. 


2.6 and can handle any number of cases of loading. 
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(4) EQFT (A, MAXA, NP, ID, NN) 

This subroutine forms the partition AF” of a substructure 

stiffness matrix, and is based on the algorithm shown in Fig. 2.7. 
(5) EQKBB (A, B, MAXA, NP, ID, NN, NLC) 

This subroutine forms the inter-boundary stiffness and load 
partitions [kK] and (Re) , and is based on the algorithm shown in Fig. 
209); 

(6) BSKB1 (B, A, MAXA, NP, NLC, ID, KB) 

This subroutine performs the backsubstitution process through 
the internal degrees of freedom, and is based on the lower algorithm of 
Pug. 2.9% 

(7) BKSK2 (A, B, MAXA, NPLI, NJ, MNLC, ID) 

This subroutine performs the backsubstitution upper algorithm 
Ohrid. 2.9). 

(8) MODAX1 (BB, MAXB, MHB, NN, ID, NWK) 

This subroutine checks the columns of the assembled stiffness 
matrix of a higher level unit stored in array BB for the actual first 
non zero component, changes array MAXB accordingly, and shifts the 
components of BB forward as may be necessary. 

(9) SKYPRD (MAXA, NPG, MHB, ID, NFIBN, NIBNS, MXND, ME) 

This subroutine forms the column heights array MHB for a 

higher level unit stiffness matrix taking into consideration the sky- 


lines of the lower level units. 


C.5.2 Description of the Arguments 


A : The stiffness storage vector. 


B : The load array. 
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BB 
ID 
KB 


LM 


MAXA 
MAXB 


ME 
MHB 
MHT 
MXND 


ND 


NFIBN 


NIBNS 


NLC 


NN 
NP 


: A higher level unit stiffness storage vector. 
: Number of degrees of freedom per node. 


: O, if the problem is unsubstructured, and 


1, if otherwise. 


: The degrees of freedom associated with an element or a 


lower level unit storage vector. 


: The diagonal component addressing array. 


: The diagonal component addressing array of a higher level 


unit. 


: Element number, or the local constituent unit number. 
: The column height array of a higher level unit. 
: The column height array. 


: The maximum number of inter-boundary nodes in any constituent 


unit. 


: Number of degrees of freedom associated with an element, or 


with the inter-boundary partition of a constituent unit. 


: Number of the first inter-boundary node in a constituent 


unit. 


: Number of inter-boundary nodes in a constituent unit. 
: Number of cases of loading. 
: Total number of nodes in a unit. 


: Total number of internal nodes in a unit. 
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C.6 The Formulation and Output Package 


This group of subroutines performs various tasks, and all are 
common to SISAPF and MUSAPF, except LOADID which is used only in MUSAPF. 
Subroutines MLOADS, DISPL, and STRESS are adapted from Reference 109. 
Descriptions of the subroutines can be found in the listing which follows. 
These subroutines are in alphabetical order. 
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